Changeset 779


Ignore:
Timestamp:
2007-12-22T18:04:17+01:00 (13 years ago)
Author:
rblod
Message:

Agrif improvment for vectorization, see ticket #41

Location:
trunk/AGRIF/AGRIF_FILES
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • trunk/AGRIF/AGRIF_FILES/modbc.F

    r774 r779  
    10461046      if (child % var % interpIndex  
    10471047     &        /= Agrif_Curgrid % parent % ngridstep ) then 
    1048        child%var%oldvalues2d(kindex:kindex+newsize-1,1)= 
    1049      &        child%var%oldvalues2d(kindex:kindex+newsize-1,2) 
     1048       child%var%oldvalues2d(1,kindex:kindex+newsize-1)= 
     1049     &        child%var%oldvalues2d(2,kindex:kindex+newsize-1) 
    10501050      endif 
    10511051 
     
    10531053      CASE (1) 
    10541054 
     1055!CDIR ALTCODE 
    10551056         do ir=bounds(1,1),bounds(1,2) 
    1056             child%var%oldvalues2d(kindex,2) = 
     1057            child%var%oldvalues2d(2,kindex) = 
    10571058     &           child%var%array1(ir) 
    10581059            kindex = kindex + 1 
     
    10621063 
    10631064        do jr=bounds(2,1),bounds(2,2) 
     1065!CDIR ALTCODE 
    10641066           do ir=bounds(1,1),bounds(1,2) 
    1065             child%var%oldvalues2d(kindex,2) = 
     1067            child%var%oldvalues2d(2,kindex) = 
    10661068     &           child%var%array2(ir,jr) 
    10671069            kindex = kindex + 1 
     
    10721074        do kr=bounds(3,1),bounds(3,2) 
    10731075          do jr=bounds(2,1),bounds(2,2) 
     1076!CDIR ALTCODE 
    10741077             do ir=bounds(1,1),bounds(1,2) 
    1075             child%var%oldvalues2d(kindex,2) = 
     1078            child%var%oldvalues2d(2,kindex) = 
    10761079     &           child%var%array3(ir,jr,kr) 
    10771080            kindex = kindex + 1 
     
    10841087           do kr=bounds(3,1),bounds(3,2) 
    10851088             do jr=bounds(2,1),bounds(2,2) 
     1089!CDIR ALTCODE 
    10861090               do ir=bounds(1,1),bounds(1,2) 
    1087             child%var%oldvalues2d(kindex,2) = 
     1091            child%var%oldvalues2d(2,kindex) = 
    10881092     &           child%var%array4(ir,jr,kr,lr) 
    10891093            kindex = kindex + 1 
     
    10981102            do kr=bounds(3,1),bounds(3,2)  
    10991103              do jr=bounds(2,1),bounds(2,2) 
     1104!CDIR ALTCODE 
    11001105                 do ir=bounds(1,1),bounds(1,2) 
    1101             child%var%oldvalues2d(kindex,2) = 
     1106            child%var%oldvalues2d(2,kindex) = 
    11021107     &           child%var%array5(ir,jr,kr,lr,mr) 
    11031108            kindex = kindex + 1 
     
    11141119               do kr=bounds(3,1),bounds(3,2) 
    11151120                do jr=bounds(2,1),bounds(2,2) 
     1121!CDIR ALTCODE 
    11161122                   do ir=bounds(1,1),bounds(1,2) 
    1117             child%var%oldvalues2d(kindex,2) = 
     1123            child%var%oldvalues2d(2,kindex) = 
    11181124     &           child%var%array6(ir,jr,kr,lr,mr,nr) 
    11191125            kindex = kindex + 1 
     
    11651171      CASE (1) 
    11661172 
     1173!CDIR ALTCODE 
    11671174         do ir=bounds(1,1),bounds(1,2) 
    11681175                child%var%array1(ir) = 
    1169      &           c2t*child % var % oldvalues2d(kindex,1)    
    1170      &         + c1t*child % var % oldvalues2d(kindex,2)     
     1176     &           c2t*child % var % oldvalues2d(1,kindex)    
     1177     &         + c1t*child % var % oldvalues2d(2,kindex)     
    11711178            kindex = kindex + 1 
    11721179         enddo        
     
    11751182 
    11761183        do jr=bounds(2,1),bounds(2,2) 
     1184!CDIR ALTCODE 
    11771185           do ir=bounds(1,1),bounds(1,2) 
    11781186                child%var%array2(ir,jr) = 
    1179      &           c2t*child % var % oldvalues2d(kindex,1)    
    1180      &         + c1t*child % var % oldvalues2d(kindex,2)  
     1187     &           c2t*child % var % oldvalues2d(1,kindex)    
     1188     &         + c1t*child % var % oldvalues2d(2,kindex)  
    11811189            kindex = kindex + 1 
    11821190           enddo 
     
    11861194        do kr=bounds(3,1),bounds(3,2) 
    11871195          do jr=bounds(2,1),bounds(2,2) 
     1196!CDIR ALTCODE 
    11881197             do ir=bounds(1,1),bounds(1,2) 
    11891198                child%var%array3(ir,jr,kr) = 
    1190      &           c2t*child % var % oldvalues2d(kindex,1)    
    1191      &         + c1t*child % var % oldvalues2d(kindex,2)  
     1199     &           c2t*child % var % oldvalues2d(1,kindex)    
     1200     &         + c1t*child % var % oldvalues2d(2,kindex)  
    11921201            kindex = kindex + 1 
    11931202             enddo 
     
    11991208           do kr=bounds(3,1),bounds(3,2) 
    12001209             do jr=bounds(2,1),bounds(2,2) 
     1210!CDIR ALTCODE 
    12011211               do ir=bounds(1,1),bounds(1,2) 
    12021212                child%var%array4(ir,jr,kr,lr) = 
    1203      &           c2t*child % var % oldvalues2d(kindex,1)    
    1204      &         + c1t*child % var % oldvalues2d(kindex,2)  
     1213     &           c2t*child % var % oldvalues2d(1,kindex)    
     1214     &         + c1t*child % var % oldvalues2d(2,kindex)  
    12051215            kindex = kindex + 1 
    12061216               enddo 
     
    12141224            do kr=bounds(3,1),bounds(3,2)  
    12151225              do jr=bounds(2,1),bounds(2,2) 
     1226!CDIR ALTCODE 
    12161227                 do ir=bounds(1,1),bounds(1,2) 
    12171228                child%var%array5(ir,jr,kr,lr,mr) = 
    1218      &           c2t*child % var % oldvalues2d(kindex,1)    
    1219      &         + c1t*child % var % oldvalues2d(kindex,2)  
     1229     &           c2t*child % var % oldvalues2d(1,kindex)    
     1230     &         + c1t*child % var % oldvalues2d(2,kindex)  
    12201231            kindex = kindex + 1 
    12211232                 enddo 
     
    12311242               do kr=bounds(3,1),bounds(3,2) 
    12321243                do jr=bounds(2,1),bounds(2,2) 
     1244!CDIR ALTCODE 
    12331245                   do ir=bounds(1,1),bounds(1,2) 
    12341246                child%var%array6(ir,jr,kr,lr,mr,nr) = 
    1235      &           c2t*child % var % oldvalues2d(kindex,1)    
    1236      &         + c1t*child % var % oldvalues2d(kindex,2)  
     1247     &           c2t*child % var % oldvalues2d(1,kindex)    
     1248     &         + c1t*child % var % oldvalues2d(2,kindex)  
    12371249            kindex = kindex + 1 
    12381250                   enddo 
     
    12771289      if (.NOT. associated(child % var % oldvalues2d)) then 
    12781290C 
    1279           allocate(child % var % oldvalues2d(newsize,2)) 
     1291          allocate(child % var % oldvalues2d(2,newsize)) 
    12801292 
    12811293          child % var % oldvalues2d=0. 
     
    12831295        else 
    12841296C 
    1285           if (SIZE(child % var % oldvalues2d,1) < newsize) then    
    1286 C 
    1287               allocate(tempoldvalues(SIZE(child % var % 
    1288      &                                    oldvalues2d,1),2)) 
     1297          if (SIZE(child % var % oldvalues2d,2) < newsize) then    
     1298C 
     1299              allocate(tempoldvalues(2,SIZE(child % var % 
     1300     &                                    oldvalues2d,2))) 
    12891301C 
    12901302              tempoldvalues = child % var % oldvalues2d 
     
    12921304              deallocate(child % var % oldvalues2d) 
    12931305C 
    1294               allocate(child % var % oldvalues2d(newsize,2)) 
     1306              allocate(child % var % oldvalues2d(2,newsize)) 
    12951307C             
    12961308              child%var%oldvalues2d=0. 
    12971309C 
    1298               child % var % oldvalues2d(1:SIZE(tempoldvalues,1),:) =  
     1310              child % var % oldvalues2d(:,1:SIZE(tempoldvalues,1)) =  
    12991311     &        tempoldvalues(:,:) 
    13001312C 
  • trunk/AGRIF/AGRIF_FILES/modbcfunction.F

    r774 r779  
    268268 
    269269        Allocate( 
    270      &    Agrif_Curgrid%tabvars(tabvarsindic)%var % oldvalues2D(1,2)) 
     270     &    Agrif_Curgrid%tabvars(tabvarsindic)%var % oldvalues2D(2,1)) 
    271271          Agrif_Curgrid%tabvars(tabvarsindic)%var % oldvalues2D = 0.  
    272272       ENDIF       
  • trunk/AGRIF/AGRIF_FILES/modinterp.F

    r662 r779  
    4242C       
    4343      IMPLICIT NONE 
     44      logical,  private:: precomputedone(7) = .FALSE. 
    4445C 
    4546      CONTAINS 
     
    13511352     &                  s_parent(nbdim),s_child(nbdim), 
    13521353     &                  ds_parent(nbdim),ds_child(nbdim),coeffraf) 
     1354                 
    13531355C                 
    13541356      Return 
     
    13991401      INTEGER i,j 
    14001402      INTEGER :: coeffraf 
     1403      REAL   , DIMENSION( 
     1404     &                pttab_child(nbdim):petab_child(nbdim), 
     1405     &                pttab_child(nbdim-1):petab_child(nbdim-1) 
     1406     &                ) :: tabout_trsp 
     1407      REAL, DIMENSION(indmin(nbdim):indmax(nbdim), 
     1408     &        pttab_child(nbdim-1):petab_child(nbdim-1)) :: tabtemp_trsp 
     1409 
    14011410C 
    14021411C 
     
    14051414C     Commentaire perso : nbdim vaut toujours 2 ici. 
    14061415C 
     1416      coeffraf = nint ( ds_parent(1) / ds_child(1) ) 
     1417      IF((TypeInterp(1) == Agrif_Linear) .AND. (coeffraf /= 1 ) )THEN 
     1418 
     1419!---CDIR NEXPAND 
     1420          IF(.NOT. precomputedone(1) ) call linear1Dprecompute2D( 
     1421     &          indmax(2)-indmin(2)+1,    
     1422     &          indmax(1)-indmin(1)+1,    
     1423     &          petab_child(1)-pttab_child(1)+1, 
     1424     &          s_parent(1),s_child(1),ds_parent(1),ds_child(1),1) 
     1425!---CDIR NEXPAND 
     1426           call linear1daftercompute(tabin,tabtemp, 
     1427     &          size(tabin), size(tabtemp),   
     1428     &          s_parent(1),s_child(1),ds_parent(1),ds_child(1),1) 
     1429  
     1430      ELSEIF((TypeInterp(1) == Agrif_PPM) .AND. (coeffraf /= 1 ) )THEN 
     1431!---CDIR NEXPAND 
     1432          IF(.NOT. precomputedone(1) ) call ppm1Dprecompute2D( 
     1433     &          indmax(2)-indmin(2)+1,    
     1434     &          indmax(1)-indmin(1)+1,    
     1435     &          petab_child(1)-pttab_child(1)+1, 
     1436     &          s_parent(1),s_child(1),ds_parent(1),ds_child(1),1) 
     1437!---CDIR NEXPAND 
     1438           call ppm1daftercompute(tabin,tabtemp, 
     1439     &          size(tabin), size(tabtemp),   
     1440     &          s_parent(1),s_child(1),ds_parent(1),ds_child(1),1) 
     1441 
     1442      ELSE 
     1443 
    14071444      do j = indmin(nbdim),indmax(nbdim) 
    14081445C         
     1446!---CDIR NEXPAND 
    14091447        Call Agrif_Interp_1D_recursive(TypeInterp(1), 
    14101448     &         tabin(indmin(nbdim-1):indmax(nbdim-1),j), 
     
    14161454C         
    14171455      enddo 
    1418 C         
     1456      ENDIF    
     1457 
    14191458      coeffraf = nint(ds_parent(nbdim)/ds_child(nbdim)) 
    14201459       
     1460      tabtemp_trsp = TRANSPOSE(tabtemp) 
     1461 
     1462      IF((TypeInterp(2) == Agrif_Linear) .AND. (coeffraf /= 1 ) )THEN 
     1463 
     1464!---CDIR NEXPAND 
     1465          IF(.NOT. precomputedone(2) ) call linear1Dprecompute2D( 
     1466     &          petab_child(1)-pttab_child(1)+1, 
     1467     &          indmax(2)-indmin(2)+1, 
     1468     &          petab_child(2)-pttab_child(2)+1, 
     1469     &          s_parent(2),s_child(2),ds_parent(2),ds_child(2),2) 
     1470!---CDIR NEXPAND 
     1471           call linear1daftercompute(tabtemp_trsp,tabout_trsp, 
     1472     &          size(tabtemp_trsp), size(tabout_trsp), 
     1473     &          s_parent(2),s_child(2),ds_parent(2),ds_child(2),2) 
     1474 
     1475      ELSEIF((TypeInterp(2) == Agrif_PPM) .AND. (coeffraf /= 1 ) )THEN 
     1476 
     1477!---CDIR NEXPAND 
     1478           IF(.NOT. precomputedone(1) )call ppm1Dprecompute2D( 
     1479     &          petab_child(1)-pttab_child(1)+1, 
     1480     &          indmax(2)-indmin(2)+1, 
     1481     &          petab_child(2)-pttab_child(2)+1, 
     1482     &          s_parent(2),s_child(2),ds_parent(2),ds_child(2),2) 
     1483!---CDIR NEXPAND 
     1484           call ppm1daftercompute(tabtemp_trsp,tabout_trsp, 
     1485     &          size(tabtemp_trsp), size(tabout_trsp), 
     1486     &          s_parent(2),s_child(2),ds_parent(2),ds_child(2),2) 
     1487 
     1488      ELSE 
    14211489      do i=pttab_child(nbdim-1),petab_child(nbdim-1) 
    14221490C 
     1491!---CDIR NEXPAND 
    14231492        Call Agrif_InterpBase(TypeInterp(2), 
    1424      &           tabtemp(i,indmin(nbdim):indmax(nbdim)), 
    1425      &                  tabout(i,pttab_child(nbdim):petab_child(nbdim)), 
     1493     &           tabtemp_trsp(indmin(nbdim):indmax(nbdim),i), 
     1494     &           tabout_trsp(pttab_child(nbdim):petab_child(nbdim),i), 
    14261495     &           indmin(nbdim),indmax(nbdim), 
    14271496     &           pttab_child(nbdim),petab_child(nbdim), 
    14281497     &           s_parent(nbdim),s_child(nbdim), 
    14291498     &           ds_parent(nbdim),ds_child(nbdim),coeffraf) 
     1499 
    14301500C         
    14311501      enddo 
     1502      ENDIF 
     1503       
     1504      tabout = TRANSPOSE(tabout_trsp) 
    14321505C 
    14331506      Return 
     
    14761549C 
    14771550C 
     1551      coeffraf = nint ( ds_parent(1) / ds_child(1) ) 
     1552      IF((TypeInterp(1) == Agrif_Linear) .AND. (coeffraf/=1))THEN  
     1553        Call linear1Dprecompute2D( 
     1554     &          indmax(2)-indmin(2)+1, 
     1555     &          indmax(1)-indmin(1)+1, 
     1556     &          petab_child(1)-pttab_child(1)+1, 
     1557     &          s_parent(1),s_child(1),ds_parent(1),ds_child(1),1) 
     1558      precomputedone(1) = .TRUE.  
     1559      ELSEIF((TypeInterp(1) == Agrif_PPM) .AND. (coeffraf/=1))THEN 
     1560        Call ppm1Dprecompute2D( 
     1561     &          indmax(2)-indmin(2)+1, 
     1562     &          indmax(1)-indmin(1)+1, 
     1563     &          petab_child(1)-pttab_child(1)+1, 
     1564     &          s_parent(1),s_child(1),ds_parent(1),ds_child(1),1) 
     1565      precomputedone(1) = .TRUE. 
     1566      ENDIF 
     1567 
     1568      coeffraf = nint ( ds_parent(2) / ds_child(2) ) 
     1569      IF((TypeInterp(2) == Agrif_Linear) .AND. (coeffraf/=1)) THEN 
     1570         Call linear1Dprecompute2D( 
     1571     &          petab_child(1)-pttab_child(1)+1, 
     1572     &          indmax(2)-indmin(2)+1, 
     1573     &          petab_child(2)-pttab_child(2)+1, 
     1574     &          s_parent(2),s_child(2),ds_parent(2),ds_child(2),2) 
     1575      precomputedone(2) = .TRUE.  
     1576      ELSEIF((TypeInterp(2) == Agrif_PPM) .AND. (coeffraf/=1)) THEN 
     1577         Call ppm1Dprecompute2D( 
     1578     &          petab_child(1)-pttab_child(1)+1, 
     1579     &          indmax(2)-indmin(2)+1, 
     1580     &          petab_child(2)-pttab_child(2)+1, 
     1581     &          s_parent(2),s_child(2),ds_parent(2),ds_child(2),2) 
     1582      precomputedone(2) = .TRUE. 
     1583      ENDIF 
     1584 
    14781585      do k = indmin(nbdim),indmax(nbdim) 
    14791586C         
     
    14901597      enddo 
    14911598       
     1599      precomputedone(1) = .FALSE. 
     1600      precomputedone(2) = .FALSE. 
     1601      coeffraf = nint ( ds_parent(3) / ds_child(3) ) 
    14921602 
    14931603      Call Agrif_Compute_nbdim_interp(s_parent(nbdim),s_child(nbdim), 
     
    18511961       ELSEIF (TYPEinterp .EQ. AGRIF_LINEAR) then     
    18521962C 
    1853 C         Linear interpolation          
     1963C         Linear interpolation  
     1964   
    18541965          Call linear1D 
    18551966     &         (parenttab,childtab, 
     
    18571968     &          s_parent,s_child,ds_parent,ds_child) 
    18581969C           
     1970      elseif ( TYPEinterp .EQ. AGRIF_PPM ) then 
     1971 
     1972          Call ppm1D 
     1973     &         (parenttab,childtab, 
     1974     &         indmax-indmin+1,petab_child-pttab_child+1, 
     1975     &         s_parent,s_child,ds_parent,ds_child) 
     1976C 
     1977 
    18591978        elseif (TYPEinterp .EQ. AGRIF_LAGRANGE) then 
    18601979C           
     
    19062025     &         s_parent,s_child,ds_parent,ds_child) 
    19072026C               
    1908       elseif ( TYPEinterp .EQ. AGRIF_PPM ) then 
    1909           Call ppm1D          
    1910      &         (parenttab,childtab, 
    1911      &         indmax-indmin+1,petab_child-pttab_child+1, 
    1912      &         s_parent,s_child,ds_parent,ds_child) 
    1913 C 
    19142027      endif  
    19152028C 
  • trunk/AGRIF/AGRIF_FILES/modinterpbasic.F

    r662 r779  
    3737C              
    3838      Real,Dimension(Agrif_MaxRaff) :: tabdiff2, tabdiff3 
    39       Real,Dimension(:),Allocatable::tabtest4       
     39      Real,Dimension(5,Agrif_MaxRaff,3) :: tabppm 
     40      Real,Dimension(:),Allocatable::tabtest4    
     41      Integer, Dimension(:,:), Allocatable :: indparent 
     42      Integer, Dimension(:,:), Allocatable ::  
     43     &           indparentppm,indchildppm 
     44      Integer, Dimension(:), Allocatable ::  
     45     &           indparentppm_1d,indchildppm_1d 
     46      Real, Dimension(:,:),Allocatable :: coeffparent    
    4047        
    4148      CONTAINS 
     
    140147C        
    141148      End Subroutine Linear1d    
     149       
     150      Subroutine Linear1dprecompute2d(np2, np,nc, 
     151     &                    s_parent,s_child,ds_parent,ds_child,dir) 
     152C 
     153CCC   Description: 
     154CCC   Subroutine to compute 2D coefficient and index for a linear 1D interpolation 
     155CCC   on a child grid (vector y) from its parent grid (vector x). 
     156C 
     157CC    Method: 
     158C 
     159C     Declarations: 
     160C 
     161 
     162C         
     163C     Arguments 
     164      INTEGER             :: np,nc,np2,dir 
     165      REAL                :: s_parent,s_child,ds_parent,ds_child 
     166C 
     167C     Local scalars 
     168      INTEGER :: i,j,coeffraf,locind_parent_left,inc,inc1,inc2,toto  
     169      Integer, Dimension(:,:), Allocatable :: indparent_tmp 
     170      Real, Dimension(:,:), Allocatable :: coeffparent_tmp 
     171      REAL    :: ypos,globind_parent_left,globind_parent_right 
     172      REAL    :: invds, invds2 
     173      REAL :: ypos2,diff 
     174C 
     175C 
     176 
     177      coeffraf = nint(ds_parent/ds_child) 
     178C 
     179      ypos = s_child 
     180 
     181      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent) 
     182 
     183        globind_parent_left = s_parent 
     184     &                        + (locind_parent_left - 1)*ds_parent 
     185 
     186        globind_parent_right = globind_parent_left + ds_parent 
     187 
     188C 
     189      invds = 1./ds_parent 
     190      invds2 = ds_child/ds_parent 
     191 
     192      ypos2 = ypos*invds 
     193      globind_parent_right=globind_parent_right*invds 
     194 
     195      if (.not.allocated(indparent)) then 
     196      allocate(indparent(nc*np2,3),coeffparent(nc*np2,3)) 
     197      else 
     198      if (size(indparent,1)<nc*np2) then 
     199      allocate(coeffparent_tmp(size(indparent,1),size(indparent,2))) 
     200      allocate(indparent_tmp(size(indparent,1),size(indparent,2))) 
     201      coeffparent_tmp=coeffparent  
     202      indparent_tmp=indparent  
     203      deallocate(indparent,coeffparent) 
     204      allocate(indparent(nc*np2,3),coeffparent(nc*np2,3)) 
     205      coeffparent(1:size(coeffparent_tmp,1),1:size(coeffparent_tmp,2)) 
     206     &  = coeffparent_tmp  
     207      indparent(1:size(indparent_tmp,1),1:size(indparent_tmp,2)) 
     208     & = indparent_tmp 
     209      deallocate(indparent_tmp,coeffparent_tmp) 
     210      endif 
     211      endif 
     212 
     213      do i = 1,nc-1 
     214C 
     215        if (ypos2 > globind_parent_right) then 
     216           locind_parent_left = locind_parent_left + 1 
     217           globind_parent_right = globind_parent_right + 1. 
     218        endif 
     219         
     220        diff=(globind_parent_right - ypos2) 
     221        indparent(i,dir) = locind_parent_left 
     222        coeffparent(i,dir) = diff 
     223         
     224        ypos2 = ypos2 + invds2 
     225C 
     226      enddo 
     227C 
     228      ypos = s_child + (nc-1)*ds_child 
     229      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent) 
     230 
     231      if (locind_parent_left == np) then 
     232      indparent(nc,dir) = locind_parent_left-1 
     233      coeffparent(nc,dir) = 0. 
     234        else 
     235C 
     236          globind_parent_left = s_parent 
     237     &                        + (locind_parent_left - 1)*ds_parent 
     238C       
     239       indparent(nc,dir) = locind_parent_left 
     240        
     241         coeffparent(nc,dir) = (globind_parent_left + ds_parent - ypos) 
     242     &       * invds 
     243      endif                                            
     244 
     245      do i=2, np2 
     246      inc  =  i*nc 
     247      inc1 = (i-1)*nc 
     248      inc2 = (i-2)*nc  
     249!CDIR ALTCODE 
     250      indparent(1+inc1:inc,dir) = indparent(1+inc2:inc1,dir)+np 
     251!CDIR ALTCODE 
     252      coeffparent(1+inc1:inc,dir) =coeffparent(1:nc,dir) 
     253      enddo       
     254 
     255      Return 
     256C 
     257C        
     258      End Subroutine Linear1dprecompute2d 
     259 
     260 
     261 
     262      Subroutine Linear1dprecompute(np,nc, 
     263     &                    s_parent,s_child,ds_parent,ds_child)  
     264C 
     265CCC   Description: 
     266CCC   Subroutine to compute 1D coefficient and index for a linear 1D interpolation 
     267CCC   on a child grid (vector y) from its parent grid (vector x). 
     268C 
     269CC    Method: 
     270C 
     271C     Declarations: 
     272C 
     273       
     274C         
     275C     Arguments 
     276      INTEGER             :: np,nc  
     277      REAL                :: s_parent,s_child,ds_parent,ds_child 
     278C 
     279C     Local scalars 
     280      INTEGER :: i,coeffraf,locind_parent_left 
     281      REAL    :: ypos,globind_parent_left,globind_parent_right 
     282      REAL    :: invds, invds2 
     283      REAL :: ypos2,diff 
     284C 
     285C 
     286 
     287      coeffraf = nint(ds_parent/ds_child) 
     288C 
     289      if (coeffraf == 1) then 
     290C 
     291          return 
     292C 
     293      endif                           
     294C 
     295      ypos = s_child  
     296 
     297      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent) 
     298 
     299        globind_parent_left = s_parent  
     300     &                        + (locind_parent_left - 1)*ds_parent 
     301 
     302        globind_parent_right = globind_parent_left + ds_parent 
     303 
     304C 
     305      invds = 1./ds_parent 
     306      invds2 = ds_child/ds_parent 
     307       
     308      ypos2 = ypos*invds 
     309      globind_parent_right=globind_parent_right*invds 
     310       
     311      if (.not.allocated(indparent)) then 
     312      allocate(indparent(nc,1),coeffparent(nc,1)) 
     313      else 
     314      if (size(indparent)<nc) then 
     315      deallocate(indparent,coeffparent) 
     316      allocate(indparent(nc,1),coeffparent(nc,1)) 
     317      endif 
     318      endif 
     319       
     320      do i = 1,nc-1 
     321C 
     322        if (ypos2 > globind_parent_right) then 
     323           locind_parent_left = locind_parent_left + 1 
     324           globind_parent_right = globind_parent_right + 1. 
     325        endif 
     326         
     327        diff=(globind_parent_right - ypos2) 
     328        indparent(i,1) = locind_parent_left 
     329        coeffparent(i,1) = diff 
     330        ypos2 = ypos2 + invds2 
     331C 
     332      enddo 
     333C 
     334      ypos = s_child + (nc-1)*ds_child 
     335      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent) 
     336       
     337      if (locind_parent_left == np) then 
     338      indparent(nc,1) = locind_parent_left-1 
     339      coeffparent(nc,1) = 0. 
     340        else 
     341C 
     342          globind_parent_left = s_parent  
     343     &                        + (locind_parent_left - 1)*ds_parent   
     344C       
     345       indparent(nc,1) = locind_parent_left 
     346        
     347         coeffparent(nc,1) = (globind_parent_left + ds_parent - ypos) 
     348     &       * invds 
     349      endif                                           
     350C            
     351      Return 
     352C 
     353C        
     354      End Subroutine Linear1dprecompute  
     355       
     356      Subroutine Linear1daftercompute(x,y,np,nc, 
     357     &                    s_parent,s_child,ds_parent,ds_child,dir)  
     358C 
     359CCC   Description: 
     360CCC   Subroutine to do a linear 1D interpolation on a child grid (vector y) from 
     361CCC   its parent grid (vector x) using precomputed coefficient and index.   
     362C 
     363CC    Method: 
     364C 
     365C     Declarations: 
     366C 
     367       
     368C         
     369C     Arguments 
     370      INTEGER             :: np,nc,dir       
     371      REAL,INTENT(IN), DIMENSION(np) :: x       
     372      REAL,INTENT(OUT), DIMENSION(nc) :: y   
     373      REAL                :: s_parent,s_child,ds_parent,ds_child 
     374C 
     375C     Local scalars 
     376      INTEGER :: i,coeffraf,locind_parent_left 
     377      REAL    :: ypos,globind_parent_left,globind_parent_right 
     378      REAL    :: invds, invds2 
     379      REAL :: ypos2,diff 
     380C 
     381C 
     382 
     383!CDIR ALTCODE 
     384!CDIR NODEP 
     385      do i = 1,nc 
     386C 
     387       y(i)=coeffparent(i,dir)*x(indparent(i,dir))+ 
     388     &      (1.-coeffparent(i,dir))*x(indparent(i,dir)+1) 
     389C 
     390      enddo                                      
     391C            
     392      Return 
     393C 
     394C        
     395      End Subroutine Linear1daftercompute             
    142396        
    143397C 
     
    569823      If (coeffraf == 1) Then 
    570824          locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent) 
     825!CDIR ALTCODE 
     826!CDIR SHORTLOOP 
    571827          y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1) 
    572828          return 
     
    582838         Allocate(tabtest4(-2*coeffraf:nc+2*coeffraf)) 
    583839         ENDIF 
    584       ENDIF   
     840      ENDIF 
     841             
    585842      ypos = s_child   
    586843C 
     
    595852C 
    596853       
     854!CDIR NOVECTOR 
    597855      Do i=1,coeffraf 
    598856        tabdiff2(i)=(real(i)-0.5)*invcoeffraf 
     
    602860      tabdiff3(1) = (1./3.)*a 
    603861      a=2.*a 
     862!CDIR NOVECTOR 
    604863      Do i=2,coeffraf 
    605864        tabdiff3(i) = tabdiff3(i-1)+(real(i)-1)*a 
     
    621880C  
    622881C 
     882!CDIR ALTCODE 
     883!CDIR SHORTLOOP 
    623884      Do i = nmin,nmax 
    624885         slope(i) = x(i) - x(i-1) 
    625886      Enddo 
    626887 
     888!CDIR ALTCODE 
     889!CDIR SHORTLOOP 
    627890      Do i = nmin+1,nmax-1 
    628891         xl(i)= 0.5*(x(i-1)+x(i)) 
     
    631894C 
    632895C apply parabolic monotonicity 
     896!CDIR ALTCODE 
     897!CDIR SHORTLOOP 
    633898       Do i = locind_parent_left,locind_parent_last 
    634899          delta(i) = xl(i+1) - xl(i) 
     
    644909        Do iparent = locind_parent_left,locind_parent_last        
    645910             pos=1 
     911!CDIR ALTCODE 
     912!CDIR SHORTLOOP 
    646913             Do jj = ipos - coeffraf/2+diffmod,ipos + coeffraf/2 
    647914C 
     
    653920             ipos = ipos + coeffraf 
    654921C 
    655         End do      
    656 C 
    657 C 
     922        End do 
     923C 
     924C 
     925!CDIR ALTCODE 
     926!CDIR SHORTLOOP 
    658927        y(1:nc)=tabtest4(1:nc) 
    659928         
    660929      Return 
    661930      End Subroutine ppm1d 
     931       
     932 
     933      Subroutine ppm1dprecompute2d(np2,np,nc, 
     934     &                    s_parent,s_child,ds_parent,ds_child,dir) 
     935C 
     936CCC   Description: 
     937CCC   Subroutine to compute 2D coefficients and index for a 1D interpolation 
     938CCC   using piecewise parabolic method   
     939CC    Method:  
     940C 
     941C     Declarations: 
     942C 
     943      Implicit none 
     944C         
     945C     Arguments 
     946      Integer             :: np2,np,nc,dir 
     947C      Real, Dimension(:),Allocatable :: ytemp 
     948      Real                :: s_parent,s_child,ds_parent,ds_child 
     949C 
     950C     Local scalars 
     951      Integer, Dimension(:,:), Allocatable :: indparent_tmp 
     952      Integer, Dimension(:,:), Allocatable :: indchild_tmp 
     953      Integer :: i,coeffraf,locind_parent_left,locind_parent_last 
     954      Integer :: iparent,ipos,pos,nmin,nmax 
     955      Real    :: ypos 
     956      integer :: i1,jj 
     957      Real :: xpmin,a 
     958C       
     959      Real :: xrmin,xrmax,am3,s2,s1 
     960      Real, Dimension(np) :: xl,delta,a6,slope 
     961      INTEGER :: diffmod 
     962      REAL :: invcoeffraf 
     963C       
     964      coeffraf = nint(ds_parent/ds_child) 
     965C 
     966      invcoeffraf = ds_child/ds_parent 
     967C  
     968 
     969      if (.not.allocated(indparentppm)) then 
     970      allocate( 
     971     & indparentppm(np2*nc,3), 
     972     & indchildppm(np2*nc,3) 
     973     &        ) 
     974      else 
     975      if (size(indparentppm,1)<np2*nc) then 
     976      allocate( 
     977     &  indparent_tmp(size(indparentppm,1),size(indparentppm,2)), 
     978     &  indchild_tmp(size(indparentppm,1),size(indparentppm,2))) 
     979      indparent_tmp = indparentppm  
     980      indchild_tmp  = indchildppm  
     981      deallocate(indparentppm,indchildppm) 
     982      allocate( 
     983     &indparentppm(np2*nc,3), 
     984     &indchildppm(np2*nc,3) 
     985     &        ) 
     986      indparentppm(1:size(indparent_tmp,1),1:size(indparent_tmp,2)) 
     987     & = indparent_tmp 
     988      indchildppm(1:size(indparent_tmp,1),1:size(indparent_tmp,2)) 
     989     & = indchild_tmp 
     990      deallocate(indparent_tmp,indchild_tmp) 
     991      endif 
     992      endif 
     993 
     994      if (.not.allocated(indparentppm_1d)) then 
     995      allocate(indparentppm_1d(-2*coeffraf:nc+2*coeffraf), 
     996     &         indchildppm_1d(-2*coeffraf:nc+2*coeffraf)) 
     997      else 
     998      if (size(indparentppm_1d)<nc+4*coeffraf+1) then 
     999      deallocate(indparentppm_1d,indchildppm_1d) 
     1000      allocate(indparentppm_1d(-2*coeffraf:nc+2*coeffraf), 
     1001     &         indchildppm_1d(-2*coeffraf:nc+2*coeffraf)) 
     1002      endif 
     1003      endif 
     1004 
     1005 
     1006      ypos = s_child 
     1007C 
     1008      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent) 
     1009      locind_parent_last = 1 + 
     1010     &      agrif_ceiling((ypos +(nc - 1) 
     1011     &      *ds_child - s_parent)/ds_parent) 
     1012C 
     1013      xpmin = s_parent + (locind_parent_left-1)*ds_parent 
     1014      i1 = 1+agrif_int((xpmin-s_child)/ds_child) 
     1015C      
     1016C 
     1017 
     1018      Do i=1,coeffraf 
     1019        tabdiff2(i)=(real(i)-0.5)*invcoeffraf 
     1020      EndDo 
     1021 
     1022      a = invcoeffraf**2 
     1023      tabdiff3(1) = (1./3.)*a 
     1024      a=2.*a 
     1025!CDIR ALTCODE 
     1026      Do i=2,coeffraf 
     1027        tabdiff3(i) = tabdiff3(i-1)+(real(i)-1)*a 
     1028      EndDo 
     1029 
     1030!CDIR ALTCODE 
     1031      Do i=1,coeffraf 
     1032      tabppm(1,i,dir) = 0.08333333333333* 
     1033     &              (-1.+4*tabdiff2(i)-3*tabdiff3(i)) 
     1034      tabppm(2,i,dir) = 0.08333333333333* 
     1035     &              (7.-26.*tabdiff2(i)+18.*tabdiff3(i)) 
     1036      tabppm(3,i,dir) = 0.08333333333333* 
     1037     &              (7.+30*tabdiff2(i)-30*tabdiff3(i)) 
     1038      tabppm(4,i,dir) = 0.08333333333333* 
     1039     &              (-1.-10.*tabdiff2(i)+18.*tabdiff3(i)) 
     1040      tabppm(5,i,dir) = 0.08333333333333* 
     1041     &              (2*tabdiff2(i)-3*tabdiff3(i))                    
     1042      End Do 
     1043C  
     1044C 
     1045        diffmod = 0 
     1046       IF (mod(coeffraf,2) == 0) diffmod = 1 
     1047C 
     1048        ipos = i1 
     1049C                
     1050        Do iparent = locind_parent_left,locind_parent_last 
     1051             pos=1 
     1052!CDIR ALTCODE 
     1053             Do jj = ipos - coeffraf/2+diffmod,ipos + coeffraf/2 
     1054         indparentppm_1d(jj) = iparent-2 
     1055         indchildppm_1d(jj) = pos 
     1056               pos = pos+1 
     1057             End do 
     1058             ipos = ipos + coeffraf 
     1059C 
     1060        End do 
     1061      
     1062       Do i=1,np2 
     1063 
     1064       indparentppm(1+(i-1)*nc:i*nc,dir) =  
     1065     &             indparentppm_1d(1:nc) + (i-1)*np 
     1066       indchildppm (1+(i-1)*nc:i*nc,dir) =  
     1067     &             indchildppm_1d (1:nc) + (i-1)*np 
     1068 
     1069       End do 
     1070 
     1071      Return 
     1072      End Subroutine ppm1dprecompute2d 
     1073 
     1074 
     1075      Subroutine ppm1dprecompute(np,nc, 
     1076     &                    s_parent,s_child,ds_parent,ds_child)  
     1077C 
     1078CCC   Description: 
     1079CCC   Subroutine to compute coefficient and index for  a 1D interpolation  
     1080CCC   using piecewise parabolic method   
     1081CC    Method: 
     1082C 
     1083C     Declarations: 
     1084C 
     1085      Implicit none 
     1086C         
     1087C     Arguments 
     1088      Integer             :: np,nc 
     1089C      Real, Dimension(:),Allocatable :: ytemp 
     1090      Real                :: s_parent,s_child,ds_parent,ds_child 
     1091C 
     1092C     Local scalars 
     1093      Integer :: i,coeffraf,locind_parent_left,locind_parent_last 
     1094      Integer :: iparent,ipos,pos,nmin,nmax 
     1095      Real    :: ypos 
     1096      integer :: i1,jj 
     1097      Real :: xpmin,a 
     1098C       
     1099      Real :: xrmin,xrmax,am3,s2,s1   
     1100      Real, Dimension(np) :: xl,delta,a6,slope 
     1101C      Real, Dimension(:),Allocatable  :: diff,diff2,diff3     
     1102      INTEGER :: diffmod 
     1103      REAL :: invcoeffraf 
     1104C       
     1105      coeffraf = nint(ds_parent/ds_child) 
     1106C 
     1107      If (coeffraf == 1) Then 
     1108          return 
     1109      End If 
     1110      invcoeffraf = ds_child/ds_parent 
     1111C  
     1112       
     1113      if (.not.allocated(indparentppm)) then 
     1114      allocate(indparentppm(-2*coeffraf:nc+2*coeffraf,1), 
     1115     &         indchildppm(-2*coeffraf:nc+2*coeffraf,1)) 
     1116      else 
     1117      if (size(indparentppm,1)<nc+4*coeffraf+1) then 
     1118      deallocate(indparentppm,indchildppm) 
     1119      allocate(indparentppm(-2*coeffraf:nc+2*coeffraf,1), 
     1120     &         indchildppm(-2*coeffraf:nc+2*coeffraf,1)) 
     1121      endif 
     1122      endif 
     1123             
     1124      ypos = s_child   
     1125C 
     1126      locind_parent_left = 1 + agrif_int((ypos - s_parent)/ds_parent)  
     1127      locind_parent_last = 1 + 
     1128     &      agrif_ceiling((ypos +(nc - 1)  
     1129     &      *ds_child - s_parent)/ds_parent)   
     1130C 
     1131      xpmin = s_parent + (locind_parent_left-1)*ds_parent        
     1132      i1 = 1+agrif_int((xpmin-s_child)/ds_child)         
     1133C      
     1134C 
     1135       
     1136      Do i=1,coeffraf 
     1137        tabdiff2(i)=(real(i)-0.5)*invcoeffraf 
     1138      EndDo 
     1139 
     1140      a = invcoeffraf**2  
     1141      tabdiff3(1) = (1./3.)*a 
     1142      a=2.*a 
     1143!CDIR ALTCODE 
     1144!!!CDIR SHORTLOOP 
     1145      Do i=2,coeffraf 
     1146        tabdiff3(i) = tabdiff3(i-1)+(real(i)-1)*a 
     1147      EndDo 
     1148       
     1149!CDIR ALTCODE 
     1150!!!CDIR SHORTLOOP 
     1151      Do i=1,coeffraf 
     1152      tabppm(1,i,1) = 0.08333333333333*(-1.+4*tabdiff2(i)-3*tabdiff3(i)) 
     1153      tabppm(2,i,1) = 0.08333333333333* 
     1154     &              (7.-26.*tabdiff2(i)+18.*tabdiff3(i)) 
     1155      tabppm(3,i,1) =0.08333333333333*(7.+30*tabdiff2(i)-30*tabdiff3(i)) 
     1156      tabppm(4,i,1) = 0.08333333333333* 
     1157     &              (-1.-10.*tabdiff2(i)+18.*tabdiff3(i))      
     1158      tabppm(5,i,1) = 0.08333333333333*(2*tabdiff2(i)-3*tabdiff3(i))  
     1159      End Do    
     1160C  
     1161C 
     1162        diffmod = 0 
     1163       IF (mod(coeffraf,2) == 0) diffmod = 1 
     1164C 
     1165        ipos = i1 
     1166C                
     1167        Do iparent = locind_parent_left,locind_parent_last        
     1168             pos=1 
     1169!CDIR ALTCODE 
     1170!CDIR SHORTLOOP 
     1171             Do jj = ipos - coeffraf/2+diffmod,ipos + coeffraf/2 
     1172         indparentppm(jj,1) = iparent-2 
     1173         indchildppm(jj,1) = pos 
     1174               pos = pos+1  
     1175             End do  
     1176             ipos = ipos + coeffraf 
     1177C 
     1178        End do 
     1179         
     1180      Return 
     1181      End Subroutine ppm1dprecompute 
     1182       
     1183      Subroutine ppm1daftercompute(x,y,np,nc, 
     1184     &                    s_parent,s_child,ds_parent,ds_child,dir)  
     1185C 
     1186CCC   Description: 
     1187CCC   Subroutine to do a 1D interpolation and apply monotonicity constraints 
     1188CCC   using piecewise parabolic method   
     1189CCC   on a child grid (vector y) from its parent grid (vector x). 
     1190CC    Method: Use precomputed coefficient and index 
     1191C 
     1192C     Declarations: 
     1193C 
     1194      Implicit none 
     1195C         
     1196C     Arguments 
     1197      Integer             :: np,nc,dir 
     1198      Real, INTENT(IN),Dimension(np) :: x       
     1199      Real, INTENT(OUT),Dimension(nc) :: y 
     1200C      Real, Dimension(:),Allocatable :: ytemp 
     1201      Real                :: s_parent,s_child,ds_parent,ds_child 
     1202C 
     1203C     Local scalars 
     1204      Integer :: i,coeffraf,locind_parent_left,locind_parent_last 
     1205      Integer :: iparent,ipos,pos,nmin,nmax 
     1206      Real    :: ypos 
     1207      integer :: i1,jj 
     1208      Real :: xpmin,a 
     1209C       
     1210      Real :: xrmin,xrmax,am3,s2,s1   
     1211      Real, Dimension(np) :: xl,delta,a6,slope 
     1212      INTEGER :: diffmod 
     1213C       
     1214C 
     1215      do i=1,nc 
     1216      y(i)=tabppm(1,indchildppm(i,dir),dir)*x(indparentppm(i,dir))+ 
     1217     &        tabppm(2,indchildppm(i,dir),dir)*x(indparentppm(i,dir)+1)+ 
     1218     &        tabppm(3,indchildppm(i,dir),dir)*x(indparentppm(i,dir)+2)+ 
     1219     &        tabppm(4,indchildppm(i,dir),dir)*x(indparentppm(i,dir)+3)+ 
     1220     &        tabppm(5,indchildppm(i,dir),dir)*x(indparentppm(i,dir)+4) 
     1221      enddo 
     1222         
     1223      Return 
     1224      End Subroutine ppm1daftercompute             
    6621225       
    6631226C     **************************************************************************   
  • trunk/AGRIF/AGRIF_FILES/modmask.F

    r662 r779  
    3535C 
    3636      IMPLICIT NONE 
    37       Integer, Parameter :: MaxSearch = 5 
     37      Integer, Parameter :: MaxSearch = 3 
    3838C 
    3939      CONTAINS 
     
    103103         IF (tempP%var%array3(i0,j0,k0) 
    104104     &                        == Agrif_SpecialValue) Then 
    105             Call CalculNewValTempP((/i0,j0,k0/), 
    106      &                             tempP,parent, 
    107      &                             ppbtab,ppetab, 
    108      &                             noraftab,nbdim) 
     105!------CDIR NEXPAND 
     106            Call CalculNewValTempP3D((/i0,j0,k0/), 
     107     &      tempP%var%array3(ppbtab(1),ppbtab(2),ppbtab(3)), 
     108     &      parent%var%array3(ppbtab(1),ppbtab(2),ppbtab(3)), 
     109     &                           ppbtab,ppetab, 
     110     &                           noraftab,MaxSearch,Agrif_SpecialValue) 
     111      
     112c            Call CalculNewValTempP((/i0,j0,k0/), 
     113c     &                             tempP,parent, 
     114c     &                             ppbtab,ppetab, 
     115c     &                             noraftab,nbdim) 
     116           
    109117         ENDIF 
    110118         enddo 
     
    118126         IF (tempP%var%array4(i0,j0,k0,l0)  
    119127     &                        == Agrif_SpecialValue) Then 
    120             Call CalculNewValTempP((/i0,j0,k0,l0/), 
    121      &                             tempP,parent, 
    122      &                             ppbtab,ppetab, 
    123      &                             noraftab,nbdim) 
     128            Call CalculNewValTempP4D((/i0,j0,k0,l0/), 
     129     &      tempP%var%array4(ppbtab(1),ppbtab(2),ppbtab(3),ppbtab(4)), 
     130     &      parent%var%array4(ppbtab(1),ppbtab(2),ppbtab(3),ppbtab(4)), 
     131     &                           ppbtab,ppetab, 
     132     &                           noraftab,MaxSearch,Agrif_SpecialValue) 
    124133         ENDIF 
    125134         enddo 
     
    221230      Valmax = min(Valmax,MaxSearch) 
    222231C 
     232!CDIR NOVECTOR 
    223233      imin = indic 
     234!CDIR NOVECTOR 
    224235      imax = indic 
    225236C 
     
    238249                  if (indic(iii).GT.ppbtab(iii)) then 
    239250 
     251!CDIR NOVECTOR 
    240252                     idecal = indic 
    241253                     idecal(iii) = idecal(iii)-1 
     
    286298            SELECT CASE(nbdim) 
    287299            CASE (1) 
     300!CDIR ALTCODE 
     301!CDIR SHORTLOOP 
    288302               do ii = imin(1),imax(1) 
    289303                    ValParent = parent%var%array1(ii) 
     
    296310            CASE (2) 
    297311               do jj = imin(2),imax(2) 
     312!CDIR ALTCODE 
     313!CDIR SHORTLOOP 
    298314               do ii = imin(1),imax(1) 
    299315                    ValParent = parent%var%array2(ii,jj) 
     
    308324               do kk = imin(3),imax(3) 
    309325               do jj = imin(2),imax(2) 
     326!CDIR ALTCODE 
     327!CDIR SHORTLOOP 
    310328               do ii = imin(1),imax(1) 
    311329                    ValParent = parent%var%array3(ii,jj,kk) 
     
    322340               do kk = imin(3),imax(3) 
    323341               do jj = imin(2),imax(2) 
     342!CDIR ALTCODE 
     343!CDIR SHORTLOOP 
    324344               do ii = imin(1),imax(1) 
    325345                    ValParent = parent%var%array4(ii,jj,kk,ll) 
     
    338358               do kk = imin(3),imax(3) 
    339359               do jj = imin(2),imax(2) 
     360!CDIR ALTCODE 
     361!CDIR SHORTLOOP 
    340362               do ii = imin(1),imax(1) 
    341363                    ValParent = parent%var%array5(ii,jj,kk,ll,mm) 
     
    356378               do kk = imin(3),imax(3) 
    357379               do jj = imin(2),imax(2) 
     380!CDIR ALTCODE 
     381!CDIR SHORTLOOP 
    358382               do ii = imin(1),imax(1) 
    359383                    ValParent = parent%var%array6(ii,jj,kk,ll,mm,nn) 
     
    412436C 
    413437C 
    414       End Module Agrif_Mask        
     438      End Module Agrif_Mask      
     439       
     440      Subroutine CalculNewValTempP3D(indic, 
     441     &               tempP,parent,ppbtab, 
     442     &               ppetab,noraftab,MaxSearch,Agrif_SpecialValue) 
     443C 
     444CCC   Description: 
     445CCC   Subroutine called in the procedure Agrif_InterpnD to recalculate the value  
     446CCC   of the parent grid variable when this one is equal to Agrif_SpecialValue.  
     447C 
     448C     Declarations: 
     449C 
     450        
     451C 
     452C     Arrays arguments       
     453      INTEGER, PARAMETER :: nbdim = 3 
     454      INTEGER,DIMENSION(nbdim) :: indic 
     455      LOGICAL,DIMENSION(nbdim) :: noraftab 
     456      INTEGER,DIMENSION(nbdim) :: ppbtab,ppetab 
     457      REAL :: Agrif_SpecialValue 
     458C 
     459C     Pointer argument 
     460      REAL,DIMENSION(ppbtab(1):ppetab(1),ppbtab(2):ppetab(2), 
     461     &    ppbtab(3):ppetab(3)) :: tempP, parent  ! Part of the parent grid used for  
     462                                      ! the interpolation of the child grid                      
     463C 
     464C     Local scalar 
     465      INTEGER                  :: i,ii,iii,jj,kk,ll,mm,nn  
     466      INTEGER,DIMENSION(nbdim) :: imin,imax,idecal 
     467      INTEGER                  :: Nbvals 
     468      REAL                     :: Res 
     469      REAL                     :: ValParent 
     470      INTEGER                  :: ValMax 
     471      INTEGER :: MaxSearch 
     472      LOGICAL :: Existunmasked 
     473C       
     474C     Local arrays       
     475C 
     476      ValMax = 1 
     477!CDIR NOVECTOR 
     478      do iii = 1 , nbdim 
     479         IF (.NOT.noraftab(iii)) THEN 
     480         ValMax = max(ValMax,ppetab(iii)-indic(iii)) 
     481         ValMax = max(ValMax,indic(iii)-ppbtab(iii)) 
     482         ENDIF 
     483      enddo 
     484C 
     485      Valmax = min(Valmax,MaxSearch) 
     486C 
     487!CDIR NOVECTOR 
     488      imin = indic 
     489!CDIR NOVECTOR 
     490      imax = indic 
     491       
     492!CDIR NOVECTOR 
     493         idecal = indic   
     494C 
     495         i = Valmax 
     496 
     497            do iii = 1 , nbdim 
     498               if (.NOT.noraftab(iii)) then 
     499                  imin(iii) = max(indic(iii) - i,ppbtab(iii)) 
     500                  imax(iii) = min(indic(iii) + i,ppetab(iii)) 
     501          
     502                  if (indic(iii).GT.ppbtab(iii)) then 
     503 
     504                     idecal(iii) = idecal(iii)-1 
     505 
     506                        if (tempP(idecal(1),idecal(2),idecal(3))  
     507     &                               == Agrif_SpecialValue) then 
     508                           imin(iii) = imax(iii) 
     509                        endif 
     510 
     511                     idecal(iii) = idecal(iii)+1 
     512                  endif 
     513 
     514               endif             
     515            enddo 
     516C 
     517            Existunmasked = .FALSE. 
     518C 
     519               do kk = imin(3),imax(3) 
     520               do jj = imin(2),imax(2) 
     521!CDIR NOVECTOR 
     522               do ii = imin(1),imax(1) 
     523                    if ( parent(ii,jj,kk) .NE. Agrif_SpecialValue) then 
     524                        Existunmasked = .TRUE. 
     525                       exit 
     526                    endif  
     527                        enddo 
     528                  enddo   
     529               enddo 
     530C 
     531C 
     532          if (.Not.Existunmasked) return 
     533C 
     534         i = 1 
     535C 
     536         do While(i <= ValMax) 
     537C 
     538 
     539            do iii = 1 , nbdim 
     540               if (.NOT.noraftab(iii)) then 
     541                  imin(iii) = max(indic(iii) - i,ppbtab(iii)) 
     542                  imax(iii) = min(indic(iii) + i,ppetab(iii)) 
     543               endif             
     544            enddo 
     545C 
     546            Res = 0. 
     547            Nbvals = 0 
     548C 
     549               do kk = imin(3),imax(3) 
     550               do jj = imin(2),imax(2) 
     551!CDIR NOVECTOR 
     552               do ii = imin(1),imax(1) 
     553                    ValParent = parent(ii,jj,kk) 
     554                    if ( ValParent .NE. Agrif_SpecialValue) then 
     555                        Res = Res + ValParent 
     556                        Nbvals = Nbvals + 1 
     557                    endif  
     558                        enddo 
     559                  enddo   
     560               enddo 
     561C 
     562C 
     563             
     564            if (Nbvals.GT.0) then 
     565              tempP(indic(1),indic(2),indic(3)) = Res/Nbvals 
     566              exit 
     567            else 
     568               i = i + 1                       
     569            endif 
     570          enddo             
     571C      
     572      End Subroutine CalculNewValTempP3D         
     573 
     574      Subroutine CalculNewValTempP4D(indic, 
     575     &               tempP,parent,ppbtab, 
     576     &               ppetab,noraftab,MaxSearch,Agrif_SpecialValue) 
     577C 
     578CCC   Description: 
     579CCC   Subroutine called in the procedure Agrif_InterpnD to recalculate the value  
     580CCC   of the parent grid variable when this one is equal to Agrif_SpecialValue.  
     581C 
     582C     Declarations: 
     583C 
     584        
     585C 
     586C     Arrays arguments       
     587      INTEGER, PARAMETER :: nbdim = 4 
     588      INTEGER,DIMENSION(nbdim) :: indic 
     589      LOGICAL,DIMENSION(nbdim) :: noraftab 
     590      INTEGER,DIMENSION(nbdim) :: ppbtab,ppetab 
     591      INTEGER :: MaxSearch       
     592      REAL :: Agrif_SpecialValue       
     593C 
     594C     Pointer argument 
     595      REAL,DIMENSION(ppbtab(1):ppetab(1),ppbtab(2):ppetab(2), 
     596     &    ppbtab(3):ppetab(3), 
     597     &    ppbtab(4):ppetab(4)) :: tempP, parent  ! Part of the parent grid used for  
     598                                      ! the interpolation of the child grid                      
     599C 
     600C     Local scalar 
     601      INTEGER                  :: i,ii,iii,jj,kk,ll,mm,nn  
     602      INTEGER,DIMENSION(nbdim) :: imin,imax,idecal 
     603      INTEGER                  :: Nbvals 
     604      REAL                     :: Res 
     605      REAL                     :: ValParent 
     606      INTEGER                  :: ValMax 
     607      LOGICAL                  :: firsttest 
     608C       
     609C     Local arrays       
     610C 
     611      ValMax = 1 
     612      do iii = 1 , nbdim 
     613         IF (.NOT.noraftab(iii)) THEN 
     614         ValMax = max(ValMax,ppetab(iii)-indic(iii)) 
     615         ValMax = max(ValMax,indic(iii)-ppbtab(iii)) 
     616         ENDIF 
     617      enddo 
     618C 
     619      Valmax = min(Valmax,MaxSearch) 
     620C 
     621      imin = indic 
     622      imax = indic 
     623C 
     624         i = 1 
     625         firsttest = .TRUE. 
     626         idecal = indic   
     627     
     628C 
     629         do While(i <= ValMax) 
     630C 
     631         IF ((i == 1).AND.(firsttest)) i = Valmax 
     632 
     633            do iii = 1 , nbdim 
     634               if (.NOT.noraftab(iii)) then 
     635                  imin(iii) = max(indic(iii) - i,ppbtab(iii)) 
     636                  imax(iii) = min(indic(iii) + i,ppetab(iii)) 
     637                  if (firsttest) then                
     638                  if (indic(iii).GT.ppbtab(iii)) then 
     639 
     640 
     641                     idecal(iii) = idecal(iii)-1 
     642            
     643                     if (tempP(idecal(1),idecal(2),idecal(3),idecal(4)) 
     644     &                               == Agrif_SpecialValue) then 
     645                           imin(iii) = imax(iii) 
     646                     endif 
     647            
     648           idecal(iii) = idecal(iii)+1 
     649                  endif 
     650                  endif 
     651               endif             
     652            enddo 
     653C 
     654            Res = 0. 
     655            Nbvals = 0 
     656C 
     657               do ll = imin(4),imax(4) 
     658               do kk = imin(3),imax(3) 
     659               do jj = imin(2),imax(2) 
     660               do ii = imin(1),imax(1) 
     661                    ValParent = parent(ii,jj,kk,ll) 
     662                    if ( ValParent .NE. Agrif_SpecialValue) then 
     663                        Res = Res + ValParent 
     664                        Nbvals = Nbvals + 1 
     665                    endif  
     666                              enddo 
     667                        enddo 
     668                  enddo   
     669               enddo 
     670C 
     671C 
     672             
     673            if (Nbvals.GT.0) then 
     674              if (firsttest) then 
     675                   firsttest = .FALSE. 
     676                   i=1 
     677                   cycle 
     678              endif 
     679          
     680              tempP(indic(1),indic(2),indic(3),indic(4)) = Res/Nbvals 
     681        
     682              exit 
     683            else 
     684               if (firsttest) exit 
     685               i = i + 1                       
     686            endif 
     687          enddo             
     688C      
     689      End Subroutine CalculNewValTempP4D         
  • trunk/AGRIF/AGRIF_FILES/modmpp.F

    r774 r779  
    146146        sendtoproc(k) = .true. 
    147147C 
     148!CDIR SHORTLOOP 
    148149        do i = 1,nbdim 
    149150C 
     
    228229                datasize = 1 
    229230C 
     231!CDIR SHORTLOOP 
    230232                do i = 1,nbdim 
    231233C 
     
    305307                datasize = 1 
    306308C 
     309!CDIR SHORTLOOP 
    307310                do i = 1,nbdim 
    308311C 
  • trunk/AGRIF/AGRIF_FILES/modupdate.F

    r662 r779  
    4242C 
    4343      IMPLICIT NONE 
     44      logical, private :: precomputedone(7) = .FALSE.       
    4445C       
    4546      CONTAINS 
     
    13011302 
    13021303      if ( nbdim .EQ. 1 ) then 
     1304         tempP%var%array1 = 0. 
    13031305         Call Agrif_Update_1D_recursive(TypeUpdate, 
    13041306     &           tempP%var%array1,tempCextend%var%array1, 
     
    18221824      REAL, DIMENSION(indmin(1):indmax(1), 
    18231825     &                 pttab_child(2):petab_child(2)) :: tabtemp 
     1826      REAL, DIMENSION(indmin(2):indmax(2), 
     1827     &                indmin(1):indmax(1))           :: tempP_trsp 
     1828      REAL, DIMENSION(pttab_child(2):petab_child(2), 
     1829     &                 indmin(1):indmax(1)) :: tabtemp_trsp 
    18241830      INTEGER :: i,j 
    18251831      INTEGER :: coeffraf,locind_child_left 
    18261832C 
     1833      tabtemp = 0. 
     1834 
     1835 
     1836      coeffraf = nint ( ds_parent(1) / ds_child(1) )  
     1837      IF((TypeUpdate(1) == AGRIF_Update_average)  
     1838     &                 .AND. (coeffraf /= 1 ))THEN 
     1839!---CDIR NEXPAND 
     1840         IF(.NOT. precomputedone(1) ) Call average1Dprecompute2D 
     1841     &       (petab_child(2)-pttab_child(2)+1, 
     1842     &        indmax(1)-indmin(1)+1, 
     1843     &        petab_child(1)-pttab_child(1)+1, 
     1844     &        s_parent(1),s_child(1),ds_parent(1),ds_child(1),1) 
     1845!---CDIR NEXPAND 
     1846          Call average1Daftercompute 
     1847     &       (  tabtemp, tempC, 
     1848     &          size(tabtemp), size(tempC), 
     1849     &          s_parent(1),s_child(1),ds_parent(1),ds_child(1),1) 
     1850  
     1851      ELSEIF((TypeUpdate(1) == AGRIF_Update_copy) 
     1852     &                 .AND. (coeffraf /= 1 ))THEN 
     1853!---CDIR NEXPAND 
     1854         IF(.NOT. precomputedone(1) ) Call copy1Dprecompute2D 
     1855     &       (petab_child(2)-pttab_child(2)+1, 
     1856     &        indmax(1)-indmin(1)+1, 
     1857     &        petab_child(1)-pttab_child(1)+1, 
     1858     &        s_parent(1),s_child(1),ds_parent(1),ds_child(1),1) 
     1859!---CDIR NEXPAND 
     1860          Call copy1Daftercompute 
     1861     &       (  tabtemp, tempC, 
     1862     &          size(tabtemp), size(tempC), 
     1863     &          s_parent(1),s_child(1),ds_parent(1),ds_child(1),1) 
     1864 
     1865      ELSE 
    18271866      do j = pttab_child(nbdim),petab_child(nbdim) 
    18281867C 
     1868!---CDIR NEXPAND 
    18291869        Call Agrif_Update_1D_recursive(TypeUpdate(1:nbdim-1),     
    18301870     &         tabtemp(:,j), 
     
    18361876C 
    18371877      enddo 
    1838        
     1878      ENDIF 
     1879      tabtemp_trsp = TRANSPOSE(tabtemp)  
     1880      coeffraf = nint(ds_parent(nbdim)/ds_child(nbdim)) 
     1881     
     1882!---CDIR NEXPAND 
    18391883      Call Agrif_Compute_nbdim_update(s_parent(nbdim),s_child(nbdim), 
    18401884     &  ds_parent(nbdim),ds_child(nbdim),coeffraf,locind_child_left) 
    18411885C 
     1886       
     1887      tempP_trsp = 0.  
     1888 
     1889      IF((TypeUpdate(2) == AGRIF_Update_average) 
     1890     &                .AND. (coeffraf /= 1 ))THEN 
     1891!---CDIR NEXPAND 
     1892         IF(.NOT. precomputedone(2) ) Call average1Dprecompute2D 
     1893     &      ( indmax(1)-indmin(1)+1, 
     1894     &        indmax(2)-indmin(2)+1, 
     1895     &        petab_child(2)-pttab_child(2)+1, 
     1896     &        s_parent(2),s_child(2),ds_parent(2),ds_child(2),2) 
     1897!---CDIR NEXPAND 
     1898          Call average1Daftercompute 
     1899     &       (  tempP_trsp, tabtemp_trsp, 
     1900     &          size(tempP_trsp), size(tabtemp_trsp), 
     1901     &          s_parent(2),s_child(2),ds_parent(2),ds_child(2),2) 
     1902 
     1903      ELSEIF((TypeUpdate(2) == AGRIF_Update_copy) 
     1904     &                .AND. (coeffraf /= 1 ))THEN 
     1905!---CDIR NEXPAND 
     1906         IF(.NOT. precomputedone(2) ) Call copy1Dprecompute2D 
     1907     &      ( indmax(1)-indmin(1)+1, 
     1908     &        indmax(2)-indmin(2)+1, 
     1909     &        petab_child(2)-pttab_child(2)+1, 
     1910     &        s_parent(2),s_child(2),ds_parent(2),ds_child(2),2) 
     1911!---CDIR NEXPAND 
     1912          Call copy1Daftercompute 
     1913     &       (  tempP_trsp, tabtemp_trsp, 
     1914     &          size(tempP_trsp), size(tabtemp_trsp), 
     1915     &          s_parent(2),s_child(2),ds_parent(2),ds_child(2),2) 
     1916 
     1917      ELSE 
     1918 
    18421919      do i = indmin(1),indmax(1) 
    18431920C 
     1921!---CDIR NEXPAND 
    18441922        Call Agrif_UpdateBase(TypeUpdate(2), 
    1845      &           tempP(i,:), 
    1846      &          tabtemp(i,:), 
     1923     &           tempP_trsp(indmin(nbdim):indmax(nbdim),i), 
     1924     &          tabtemp_trsp(pttab_child(nbdim):petab_child(nbdim),i), 
    18471925     &           indmin(nbdim),indmax(nbdim), 
    18481926     &           pttab_child(nbdim),petab_child(nbdim), 
     
    18521930C         
    18531931      enddo 
     1932       
     1933      ENDIF 
     1934 
     1935      tempP = TRANSPOSE(tempP_trsp) 
    18541936C 
    18551937      Return 
     
    19021984C 
    19031985C 
     1986      coeffraf = nint ( ds_parent(1) / ds_child(1) ) 
     1987      IF((TypeUpdate(1) == AGRIF_Update_average) 
     1988     &                 .AND. (coeffraf /= 1 ))THEN 
     1989!---CDIR NEXPAND 
     1990         Call average1Dprecompute2D 
     1991     &       (petab_child(2)-pttab_child(2)+1, 
     1992     &        indmax(1)-indmin(1)+1, 
     1993     &        petab_child(1)-pttab_child(1)+1, 
     1994     &        s_parent(1),s_child(1),ds_parent(1),ds_child(1),1) 
     1995      precomputedone(1) = .TRUE. 
     1996      ELSEIF((TypeUpdate(1) == AGRIF_Update_copy) 
     1997     &                 .AND. (coeffraf /= 1 ))THEN 
     1998!---CDIR NEXPAND 
     1999         Call copy1Dprecompute2D 
     2000     &       (petab_child(2)-pttab_child(2)+1, 
     2001     &        indmax(1)-indmin(1)+1, 
     2002     &        petab_child(1)-pttab_child(1)+1, 
     2003     &        s_parent(1),s_child(1),ds_parent(1),ds_child(1),1) 
     2004      precomputedone(1) = .TRUE. 
     2005      ENDIF 
     2006 
     2007      coeffraf = nint ( ds_parent(2) / ds_child(2) ) 
     2008      IF((TypeUpdate(2) == AGRIF_Update_average) 
     2009     &                .AND. (coeffraf /= 1 ))THEN 
     2010!---CDIR NEXPAND 
     2011         Call average1Dprecompute2D 
     2012     &      ( indmax(1)-indmin(1)+1, 
     2013     &        indmax(2)-indmin(2)+1, 
     2014     &        petab_child(2)-pttab_child(2)+1, 
     2015     &        s_parent(2),s_child(2),ds_parent(2),ds_child(2),2) 
     2016      precomputedone(2) = .TRUE. 
     2017      ELSEIF((TypeUpdate(2) == AGRIF_Update_copy) 
     2018     &                .AND. (coeffraf /= 1 ))THEN 
     2019!---CDIR NEXPAND 
     2020         Call copy1Dprecompute2D 
     2021     &      ( indmax(1)-indmin(1)+1, 
     2022     &        indmax(2)-indmin(2)+1, 
     2023     &        petab_child(2)-pttab_child(2)+1, 
     2024     &        s_parent(2),s_child(2),ds_parent(2),ds_child(2),2) 
     2025      precomputedone(2) = .TRUE. 
     2026      ENDIF 
     2027 
     2028 
    19042029      do k = pttab_child(nbdim),petab_child(nbdim) 
    19052030C 
     
    19142039      enddo 
    19152040C 
     2041      precomputedone(1) = .FALSE. 
     2042      precomputedone(2) = .FALSE. 
     2043      coeffraf = nint ( ds_parent(3) / ds_child(3) ) 
     2044 
    19162045      Call Agrif_Compute_nbdim_update(s_parent(nbdim),s_child(nbdim), 
    19172046     &  ds_parent(nbdim),ds_child(nbdim),coeffraf,locind_child_left) 
     
    19302059               
    19312060      ELSE 
     2061      tempP = 0. 
    19322062C 
    19332063      do j = indmin(2),indmax(2) 
     
    20212151     &  ds_parent(nbdim),ds_child(nbdim),coeffraf,locind_child_left) 
    20222152C 
     2153      tempP = 0. 
     2154 
    20232155      do k = indmin(3),indmax(3) 
    20242156C 
     
    21202252      Call Agrif_Compute_nbdim_update(s_parent(nbdim),s_child(nbdim), 
    21212253     &  ds_parent(nbdim),ds_child(nbdim),coeffraf,locind_child_left) 
     2254      tempP = 0. 
    21222255C 
    21232256      do l = indmin(4),indmax(4) 
     
    22322365     &  ds_parent(nbdim),ds_child(nbdim),coeffraf,locind_child_left) 
    22332366C 
     2367      tempP = 0. 
     2368 
    22342369      do m = indmin(5),indmax(5) 
    22352370      do l = indmin(4),indmax(4) 
     
    23072442C 
    23082443        elseif (TypeUpdate == AGRIF_Update_average) then 
    2309 C              
     2444C            
    23102445          Call average1D 
    23112446     &       (parenttab,childtab, 
    23122447     &          indmax-indmin+1,petab_child-pttab_child+1, 
    2313      &          s_parent,s_child,ds_parent,ds_child)     
     2448     &          s_parent,s_child,ds_parent,ds_child)  
    23142449C 
    23152450        elseif (TypeUpdate == AGRIF_Update_full_weighting) then 
     
    24212556#endif 
    24222557            
    2423       End Module Agrif_Update 
    24242558 
    24252559C     ************************************************************************** 
     
    24432577C 
    24442578C     Arguments 
    2445       USE Agrif_Update 
    24462579      INTEGER                   :: nbdim 
    24472580      INTEGER, DIMENSION(nbdim) :: TypeUpdate            ! TYPE of update (copy or average) 
     
    24662599     &                  ds_parent(nbdim),ds_child(nbdim), 
    24672600     &                  coeffraf,locind_child_left) 
     2601 
    24682602C 
    24692603      Return 
     
    24712605C 
    24722606      End Subroutine Agrif_Update_1D_recursive       
     2607 
     2608      End Module Agrif_Update 
  • trunk/AGRIF/AGRIF_FILES/modupdatebasic.F

    r662 r779  
    3737       
    3838      IMPLICIT NONE 
     39       
     40      integer,dimension(:,:),allocatable :: indchildcopy 
     41      integer,dimension(:,:),allocatable :: indchildaverage 
    3942C              
    4043 
     
    7881          locind_child_left = 1 + nint((s_parent - s_child)/ds_child) 
    7982C         
     83!CDIR ALTCODE 
    8084          x(1:np) = y(locind_child_left:locind_child_left+np-1) 
    8185C 
     
    8791      locind_child_left = 1 + nint((s_parent - s_child)/ds_child) 
    8892       
     93!CDIR ALTCODE 
    8994      do i = 1,np 
    9095C       
     
    100105C 
    101106      End Subroutine agrif_copy1d 
    102 C 
    103 C 
     107 
     108C     **************************************************************************   
     109CCC   Subroutine Copy1dprecompute   
     110C     **************************************************************************  
     111C 
     112      Subroutine copy1dprecompute2d(nc2,np,nc, 
     113     &                  s_parent,s_child,ds_parent,ds_child,dir) 
     114C 
     115CCC   Description: 
     116CCC   Subroutine to precompute index for a copy on a parent grid (vector x) from 
     117CCC   its child grid (vector y). 
     118C 
     119CC    Method: 
     120C 
     121C     Declarations: 
     122C 
     123 
     124C         
     125C     Arguments 
     126      INTEGER             :: nc2,np,nc 
     127      INTEGER             :: dir    
     128      REAL                :: s_parent,s_child 
     129      REAL                :: ds_parent,ds_child 
     130C 
     131C     Local variables 
     132      INTEGER,DIMENSION(:,:),ALLOCATABLE :: indchildcopy_tmp 
     133      INTEGER :: i,locind_child_left,coeffraf 
     134C  
     135C 
     136      coeffraf = nint(ds_parent/ds_child) 
     137C 
     138      locind_child_left = 1 + nint((s_parent - s_child)/ds_child) 
     139 
     140      if (.not.allocated(indchildcopy)) then 
     141      allocate(indchildcopy(np*nc2,3)) 
     142      else 
     143      if (size(indchildcopy,1)<np*nc2) then 
     144      allocate( indchildcopy_tmp( 
     145     &         size(indchildcopy,1),size(indchildcopy,2))) 
     146      indchildcopy_tmp = indchildcopy 
     147      deallocate(indchildcopy) 
     148      allocate(indchildcopy(np*nc2,3)) 
     149      indchildcopy(1:size(indchildcopy_tmp,1), 
     150     &              1:size(indchildcopy_tmp,2)) = indchildcopy_tmp 
     151      deallocate(indchildcopy_tmp) 
     152      endif 
     153      endif 
     154 
     155 
     156      do i = 1,np 
     157C       
     158         indchildcopy(i,dir) = locind_child_left 
     159         locind_child_left = locind_child_left + coeffraf 
     160C          
     161      enddo 
     162 
     163      do i =2, nc2 
     164        indchildcopy(1+(i-1)*np:i*np,dir)= 
     165     &       indchildcopy(1+(i-2)*np:(i-1)*np,dir) + nc  
     166      enddo 
     167 
     168C 
     169      Return 
     170C 
     171C 
     172      End Subroutine copy1dprecompute2d 
     173 
     174C 
     175C 
     176C     **************************************************************************   
     177CCC   Subroutine Copy1dprecompute   
     178C     **************************************************************************  
     179C 
     180      Subroutine copy1dprecompute(np,nc, 
     181     &                  s_parent,s_child,ds_parent,ds_child,dir) 
     182C 
     183CCC   Description: 
     184CCC   Subroutine to precompute index for a copy on a parent grid (vector x) from 
     185CCC   its child grid (vector y). 
     186C 
     187CC    Method: 
     188C 
     189C     Declarations: 
     190C 
     191 
     192C         
     193C     Arguments 
     194      INTEGER             :: np,nc 
     195      INTEGER             :: dir 
     196      REAL                :: s_parent,s_child 
     197      REAL                :: ds_parent,ds_child 
     198C 
     199C     Local variables 
     200      INTEGER :: i,locind_child_left,coeffraf 
     201C  
     202C 
     203      coeffraf = nint(ds_parent/ds_child) 
     204C 
     205      locind_child_left = 1 + nint((s_parent - s_child)/ds_child) 
     206 
     207      if (.not.allocated(indchildcopy)) then 
     208      allocate(indchildcopy(np,1)) 
     209      else 
     210      if (size(indchildcopy)<np) then 
     211      deallocate(indchildcopy) 
     212      allocate(indchildcopy(np,1)) 
     213      endif 
     214      endif 
     215 
     216!CDIR ALTCODE 
     217      do i = 1,np 
     218C       
     219         indchildcopy(i,1) = locind_child_left 
     220         locind_child_left = locind_child_left + coeffraf 
     221C          
     222      enddo 
     223 
     224C 
     225      Return 
     226C 
     227C 
     228      End Subroutine copy1dprecompute 
     229C 
     230C 
     231C     **************************************************************************   
     232CCC   Subroutine Copy1daftercompute   
     233C     **************************************************************************  
     234C 
     235      Subroutine copy1daftercompute(x,y,np,nc, 
     236     &                  s_parent,s_child,ds_parent,ds_child,dir) 
     237C 
     238CCC   Description: 
     239CCC   Subroutine to do a copy on a parent grid (vector x) from its child grid  
     240CCC   (vector y) using precomputed index.   
     241C 
     242CC    Method: 
     243C 
     244C     Declarations: 
     245C 
     246   
     247C         
     248C     Arguments 
     249      INTEGER             :: np,nc 
     250      INTEGER             :: dir 
     251      REAL, DIMENSION(np) :: x 
     252      REAL, DIMENSION(nc) :: y 
     253      REAL                :: s_parent,s_child 
     254      REAL                :: ds_parent,ds_child 
     255C 
     256C     Local variables 
     257      INTEGER :: i 
     258C  
     259C 
     260       
     261!CDIR ALTCODE 
     262      do i = 1,np 
     263C       
     264         x(i) = y(indchildcopy(i,dir)) 
     265C 
     266      enddo 
     267   
     268C 
     269      Return 
     270C 
     271C 
     272      End Subroutine copy1daftercompute 
     273 
     274 
    104275C 
    105276C     **************************************************************************   
     
    123294C     Local variables 
    124295      INTEGER :: i,locind_child_left,coeffraf,ii 
    125       REAL    :: xpos  
     296      REAL    :: xpos, invcoeffraf  
    126297      INTEGER :: nbnonnuls 
    127298      INTEGER :: diffmod 
     
    129300C 
    130301      coeffraf = nint(ds_parent/ds_child) 
     302      invcoeffraf = 1./coeffraf 
    131303C 
    132304      if (coeffraf == 1) then 
     
    149321 
    150322        locind_child_left = 1 + agrif_int((xpos - s_child)/ds_child) 
    151        
    152       do i = 1,np 
    153 C 
    154  
    155 C 
    156 c        if ((locind_child_left-1 < 1)  
    157 c     &      .OR. (locind_child_left+1 > nc)) then  
    158 C 
    159 c            x(i) = y(locind_child_left)                 
    160 C 
    161 c          else  
     323C 
     324      IF (Agrif_UseSpecialValueInUpdate) THEN 
     325 
     326      do i = 1,np 
    162327          nbnonnuls = 0 
     328!CDIR NOVECTOR 
    163329          Do ii = -coeffraf/2+locind_child_left+diffmod, 
    164330     &                coeffraf/2+locind_child_left 
    165331C  
    166             IF (Agrif_UseSpecialValueInUpdate) THEN 
    167             IF (y(ii) .NE. Agrif_SpecialValueFineGrid) THEN 
     332               IF (y(ii) .NE. Agrif_SpecialValueFineGrid) THEN 
    168333               nbnonnuls = nbnonnuls + 1 
    169334               x(i) = x(i) + y(ii) 
    170             ENDIF 
    171             ELSE 
     335               ENDIF 
     336          End Do 
     337               IF (nbnonnuls .NE. 0) THEN 
     338                  x(i) = x(i)/nbnonnuls 
     339               ELSE 
     340                  x(i) = Agrif_SpecialValueFineGrid 
     341               ENDIF 
     342        locind_child_left = locind_child_left + coeffraf 
     343      enddo 
     344 
     345      ELSE  
     346 
     347!CDIR ALTCODE 
     348      do i = 1,np 
     349!CDIR NOVECTOR 
     350          Do ii = -coeffraf/2+locind_child_left+diffmod, 
     351     &                coeffraf/2+locind_child_left 
     352C  
    172353               x(i) = x(i) + y(ii) 
    173             ENDIF 
    174354          End Do 
    175             IF (Agrif_UseSpecialValueInUpdate) THEN 
    176                  IF (nbnonnuls .NE. 0) THEN 
    177                     x(i) = x(i)/nbnonnuls 
    178                  ELSE 
    179                     x(i) = Agrif_SpecialValueFineGrid 
    180                  ENDIF 
    181             ELSE 
    182                  x(i) = x(i)/coeffraf 
    183             ENDIF 
    184 C 
    185 c       endif 
    186 C 
    187 c        xpos = xpos + ds_parent 
     355                 x(i) = x(i)*invcoeffraf 
    188356        locind_child_left = locind_child_left + coeffraf 
    189 C 
    190       enddo 
    191 C 
     357      enddo 
     358 
     359      ENDIF 
     360 
    192361      Return  
    193362C             
    194363C      
    195364      End Subroutine average1d 
     365       
     366      Subroutine average1dprecompute2d(nc2,np,nc, 
     367     &                     s_parent,s_child,ds_parent,ds_child,dir) 
     368C 
     369CCC   Description: 
     370CCC   Subroutine to do an update by average on a parent grid (vector x)from its  
     371CCC   child grid (vector y). 
     372C 
     373C     Arguments 
     374      INTEGER             :: nc2,np,nc,dir 
     375      REAL                :: s_parent,s_child 
     376      REAL                :: ds_parent,ds_child 
     377C 
     378C     Local variables 
     379      INTEGER :: i,locind_child_left,coeffraf,ii 
     380      INTEGER,DIMENSION(:,:),ALLOCATABLE :: indchildaverage_tmp 
     381      REAL    :: xpos, invcoeffraf 
     382      INTEGER :: nbnonnuls 
     383      INTEGER :: diffmod 
     384C  
     385C 
     386      coeffraf = nint(ds_parent/ds_child) 
     387C 
     388      xpos = s_parent 
     389C 
     390      diffmod = 0 
     391 
     392      IF ( mod(coeffraf,2) == 0 ) diffmod = 1 
     393 
     394      locind_child_left = 1 + agrif_int((xpos - s_child)/ds_child) 
     395 
     396      if (.not.allocated(indchildaverage)) then 
     397      allocate(indchildaverage(np*nc2,3)) 
     398      else 
     399      if (size(indchildaverage,1)<np*nc2) then 
     400      allocate( indchildaverage_tmp( 
     401     &         size(indchildaverage,1),size(indchildaverage,2))) 
     402      indchildaverage_tmp = indchildaverage 
     403      deallocate(indchildaverage) 
     404      allocate(indchildaverage(np*nc2,3)) 
     405      indchildaverage(1:size(indchildaverage_tmp,1), 
     406     &              1:size(indchildaverage_tmp,2)) = indchildaverage_tmp 
     407      deallocate(indchildaverage_tmp) 
     408      endif 
     409      endif 
     410 
     411      do i = 1,np 
     412        indchildaverage(i,dir)= -coeffraf/2+locind_child_left 
     413     &                                    +diffmod 
     414        locind_child_left = locind_child_left + coeffraf 
     415      enddo 
     416      do i = 2,nc2 
     417        indchildaverage(1+(i-1)*np:i*np,dir)=  
     418     &                 indchildaverage(1+(i-2)*np:(i-1)*np,dir) + nc  
     419      enddo 
     420 
     421      Return  
     422C             
     423C      
     424      End Subroutine average1dprecompute2d 
     425       
     426 
     427      Subroutine average1dprecompute(np,nc, 
     428     &                     s_parent,s_child,ds_parent,ds_child)  
     429C 
     430CCC   Description: 
     431CCC   Subroutine to do an update by average on a parent grid (vector x)from its  
     432CCC   child grid (vector y). 
     433C 
     434C     Arguments 
     435      INTEGER             :: np,nc  
     436      REAL                :: s_parent,s_child 
     437      REAL                :: ds_parent,ds_child 
     438C 
     439C     Local variables 
     440      INTEGER :: i,locind_child_left,coeffraf,ii 
     441      REAL    :: xpos, invcoeffraf  
     442      INTEGER :: nbnonnuls 
     443      INTEGER :: diffmod 
     444C  
     445C 
     446      coeffraf = nint(ds_parent/ds_child) 
     447 
     448C 
     449      if (coeffraf == 1) then 
     450C 
     451          return 
     452C 
     453      endif 
     454             
     455C 
     456      xpos = s_parent 
     457C 
     458      diffmod = 0 
     459       
     460      IF ( mod(coeffraf,2) == 0 ) diffmod = 1 
     461 
     462        locind_child_left = 1 + agrif_int((xpos - s_child)/ds_child) 
     463 
     464      if (.not.allocated(indchildaverage)) then 
     465      allocate(indchildaverage(np,1)) 
     466      else 
     467      if (size(indchildaverage,1)<np) then 
     468      deallocate(indchildaverage) 
     469      allocate(indchildaverage(np,1)) 
     470      endif 
     471      endif 
     472       
     473!CDIR ALTCODE 
     474      do i = 1,np 
     475C  
     476        indchildaverage(i,1)=-coeffraf/2+locind_child_left+diffmod 
     477        locind_child_left = locind_child_left + coeffraf 
     478      enddo 
     479 
     480      Return  
     481C             
     482C      
     483      End Subroutine average1dprecompute 
     484       
     485      Subroutine average1daftercompute(x,y,np,nc, 
     486     &                     s_parent,s_child,ds_parent,ds_child,dir)  
     487C 
     488CCC   Description: 
     489CCC   Subroutine to do an update by average on a parent grid (vector x)from its  
     490CCC   child grid (vector y). 
     491C 
     492C     Arguments 
     493      INTEGER             :: np,nc,dir       
     494      REAL, DIMENSION(np) :: x       
     495      REAL, DIMENSION(nc) :: y   
     496      REAL                :: s_parent,s_child 
     497      REAL                :: ds_parent,ds_child 
     498C 
     499C     Local variables 
     500      INTEGER :: i,locind_child_left,coeffraf,ii,j 
     501      REAL    :: xpos, invcoeffraf  
     502      REAL, PARAMETER :: one_third=1./3. 
     503      INTEGER, DIMENSION(np) :: nbnonnuls 
     504      INTEGER :: diffmod 
     505      REAL, DIMENSION(0:5) :: invcoeff=(/1.,1.,0.5,one_third,0.25,0.2/) 
     506C  
     507C 
     508      coeffraf = nint(ds_parent/ds_child) 
     509      invcoeffraf = 1./coeffraf       
     510C       
     511C 
     512      IF (Agrif_UseSpecialValueInUpdate) THEN 
     513 
     514      nbnonnuls = 0 
     515      do  j =1, coeffraf 
     516          do i=1, np 
     517               IF (y(indchildaverage(i,dir) + j -1) .NE.  
     518     &               Agrif_SpecialValueFineGrid) THEN 
     519               nbnonnuls(i) = nbnonnuls(i) + 1 
     520               x(i) = x(i) +  y(indchildaverage(i,dir) + j -1 ) 
     521               ENDIF 
     522          End Do 
     523      enddo 
     524      do i=1,np 
     525      x(i)= x(i)*invcoeff(nbnonnuls(i))  
     526      enddo 
     527 
     528      ELSE 
     529!CDIR NOLOOPCHG 
     530      do  j =1, coeffraf 
     531!CDIR VECTOR 
     532          do i=1, np 
     533             x(i) = x(i) + y(indchildaverage(i,dir) + j -1 ) 
     534          enddo 
     535      enddo 
     536       x = x * invcoeffraf 
     537      ENDIF 
     538 
     539      Return  
     540C             
     541C      
     542      End Subroutine average1daftercompute             
    196543C 
    197544C 
Note: See TracChangeset for help on using the changeset viewer.