Changeset 779
- Timestamp:
- 2007-12-22T18:04:17+01:00 (16 years ago)
- Location:
- trunk/AGRIF/AGRIF_FILES
- Files:
-
- 8 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/AGRIF/AGRIF_FILES/modbc.F
r774 r779 1046 1046 if (child % var % interpIndex 1047 1047 & /= 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) 1050 1050 endif 1051 1051 … … 1053 1053 CASE (1) 1054 1054 1055 !CDIR ALTCODE 1055 1056 do ir=bounds(1,1),bounds(1,2) 1056 child%var%oldvalues2d( kindex,2) =1057 child%var%oldvalues2d(2,kindex) = 1057 1058 & child%var%array1(ir) 1058 1059 kindex = kindex + 1 … … 1062 1063 1063 1064 do jr=bounds(2,1),bounds(2,2) 1065 !CDIR ALTCODE 1064 1066 do ir=bounds(1,1),bounds(1,2) 1065 child%var%oldvalues2d( kindex,2) =1067 child%var%oldvalues2d(2,kindex) = 1066 1068 & child%var%array2(ir,jr) 1067 1069 kindex = kindex + 1 … … 1072 1074 do kr=bounds(3,1),bounds(3,2) 1073 1075 do jr=bounds(2,1),bounds(2,2) 1076 !CDIR ALTCODE 1074 1077 do ir=bounds(1,1),bounds(1,2) 1075 child%var%oldvalues2d( kindex,2) =1078 child%var%oldvalues2d(2,kindex) = 1076 1079 & child%var%array3(ir,jr,kr) 1077 1080 kindex = kindex + 1 … … 1084 1087 do kr=bounds(3,1),bounds(3,2) 1085 1088 do jr=bounds(2,1),bounds(2,2) 1089 !CDIR ALTCODE 1086 1090 do ir=bounds(1,1),bounds(1,2) 1087 child%var%oldvalues2d( kindex,2) =1091 child%var%oldvalues2d(2,kindex) = 1088 1092 & child%var%array4(ir,jr,kr,lr) 1089 1093 kindex = kindex + 1 … … 1098 1102 do kr=bounds(3,1),bounds(3,2) 1099 1103 do jr=bounds(2,1),bounds(2,2) 1104 !CDIR ALTCODE 1100 1105 do ir=bounds(1,1),bounds(1,2) 1101 child%var%oldvalues2d( kindex,2) =1106 child%var%oldvalues2d(2,kindex) = 1102 1107 & child%var%array5(ir,jr,kr,lr,mr) 1103 1108 kindex = kindex + 1 … … 1114 1119 do kr=bounds(3,1),bounds(3,2) 1115 1120 do jr=bounds(2,1),bounds(2,2) 1121 !CDIR ALTCODE 1116 1122 do ir=bounds(1,1),bounds(1,2) 1117 child%var%oldvalues2d( kindex,2) =1123 child%var%oldvalues2d(2,kindex) = 1118 1124 & child%var%array6(ir,jr,kr,lr,mr,nr) 1119 1125 kindex = kindex + 1 … … 1165 1171 CASE (1) 1166 1172 1173 !CDIR ALTCODE 1167 1174 do ir=bounds(1,1),bounds(1,2) 1168 1175 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) 1171 1178 kindex = kindex + 1 1172 1179 enddo … … 1175 1182 1176 1183 do jr=bounds(2,1),bounds(2,2) 1184 !CDIR ALTCODE 1177 1185 do ir=bounds(1,1),bounds(1,2) 1178 1186 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) 1181 1189 kindex = kindex + 1 1182 1190 enddo … … 1186 1194 do kr=bounds(3,1),bounds(3,2) 1187 1195 do jr=bounds(2,1),bounds(2,2) 1196 !CDIR ALTCODE 1188 1197 do ir=bounds(1,1),bounds(1,2) 1189 1198 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) 1192 1201 kindex = kindex + 1 1193 1202 enddo … … 1199 1208 do kr=bounds(3,1),bounds(3,2) 1200 1209 do jr=bounds(2,1),bounds(2,2) 1210 !CDIR ALTCODE 1201 1211 do ir=bounds(1,1),bounds(1,2) 1202 1212 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) 1205 1215 kindex = kindex + 1 1206 1216 enddo … … 1214 1224 do kr=bounds(3,1),bounds(3,2) 1215 1225 do jr=bounds(2,1),bounds(2,2) 1226 !CDIR ALTCODE 1216 1227 do ir=bounds(1,1),bounds(1,2) 1217 1228 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) 1220 1231 kindex = kindex + 1 1221 1232 enddo … … 1231 1242 do kr=bounds(3,1),bounds(3,2) 1232 1243 do jr=bounds(2,1),bounds(2,2) 1244 !CDIR ALTCODE 1233 1245 do ir=bounds(1,1),bounds(1,2) 1234 1246 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) 1237 1249 kindex = kindex + 1 1238 1250 enddo … … 1277 1289 if (.NOT. associated(child % var % oldvalues2d)) then 1278 1290 C 1279 allocate(child % var % oldvalues2d( newsize,2))1291 allocate(child % var % oldvalues2d(2,newsize)) 1280 1292 C 1281 1293 child % var % oldvalues2d=0. … … 1283 1295 else 1284 1296 C 1285 if (SIZE(child % var % oldvalues2d, 1) < newsize) then1286 C 1287 allocate(tempoldvalues( SIZE(child % var %1288 & oldvalues2d, 1),2))1297 if (SIZE(child % var % oldvalues2d,2) < newsize) then 1298 C 1299 allocate(tempoldvalues(2,SIZE(child % var % 1300 & oldvalues2d,2))) 1289 1301 C 1290 1302 tempoldvalues = child % var % oldvalues2d … … 1292 1304 deallocate(child % var % oldvalues2d) 1293 1305 C 1294 allocate(child % var % oldvalues2d( newsize,2))1306 allocate(child % var % oldvalues2d(2,newsize)) 1295 1307 C 1296 1308 child%var%oldvalues2d=0. 1297 1309 C 1298 child % var % oldvalues2d( 1:SIZE(tempoldvalues,1),:) =1310 child % var % oldvalues2d(:,1:SIZE(tempoldvalues,1)) = 1299 1311 & tempoldvalues(:,:) 1300 1312 C -
trunk/AGRIF/AGRIF_FILES/modbcfunction.F
r774 r779 268 268 269 269 Allocate( 270 & Agrif_Curgrid%tabvars(tabvarsindic)%var % oldvalues2D( 1,2))270 & Agrif_Curgrid%tabvars(tabvarsindic)%var % oldvalues2D(2,1)) 271 271 Agrif_Curgrid%tabvars(tabvarsindic)%var % oldvalues2D = 0. 272 272 ENDIF -
trunk/AGRIF/AGRIF_FILES/modinterp.F
r662 r779 42 42 C 43 43 IMPLICIT NONE 44 logical, private:: precomputedone(7) = .FALSE. 44 45 C 45 46 CONTAINS … … 1351 1352 & s_parent(nbdim),s_child(nbdim), 1352 1353 & ds_parent(nbdim),ds_child(nbdim),coeffraf) 1354 1353 1355 C 1354 1356 Return … … 1399 1401 INTEGER i,j 1400 1402 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 1401 1410 C 1402 1411 C … … 1405 1414 C Commentaire perso : nbdim vaut toujours 2 ici. 1406 1415 C 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 1407 1444 do j = indmin(nbdim),indmax(nbdim) 1408 1445 C 1446 !---CDIR NEXPAND 1409 1447 Call Agrif_Interp_1D_recursive(TypeInterp(1), 1410 1448 & tabin(indmin(nbdim-1):indmax(nbdim-1),j), … … 1416 1454 C 1417 1455 enddo 1418 C 1456 ENDIF 1457 1419 1458 coeffraf = nint(ds_parent(nbdim)/ds_child(nbdim)) 1420 1459 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 1421 1489 do i=pttab_child(nbdim-1),petab_child(nbdim-1) 1422 1490 C 1491 !---CDIR NEXPAND 1423 1492 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), 1426 1495 & indmin(nbdim),indmax(nbdim), 1427 1496 & pttab_child(nbdim),petab_child(nbdim), 1428 1497 & s_parent(nbdim),s_child(nbdim), 1429 1498 & ds_parent(nbdim),ds_child(nbdim),coeffraf) 1499 1430 1500 C 1431 1501 enddo 1502 ENDIF 1503 1504 tabout = TRANSPOSE(tabout_trsp) 1432 1505 C 1433 1506 Return … … 1476 1549 C 1477 1550 C 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 1478 1585 do k = indmin(nbdim),indmax(nbdim) 1479 1586 C … … 1490 1597 enddo 1491 1598 1599 precomputedone(1) = .FALSE. 1600 precomputedone(2) = .FALSE. 1601 coeffraf = nint ( ds_parent(3) / ds_child(3) ) 1492 1602 1493 1603 Call Agrif_Compute_nbdim_interp(s_parent(nbdim),s_child(nbdim), … … 1851 1961 ELSEIF (TYPEinterp .EQ. AGRIF_LINEAR) then 1852 1962 C 1853 C Linear interpolation 1963 C Linear interpolation 1964 1854 1965 Call linear1D 1855 1966 & (parenttab,childtab, … … 1857 1968 & s_parent,s_child,ds_parent,ds_child) 1858 1969 C 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) 1976 C 1977 1859 1978 elseif (TYPEinterp .EQ. AGRIF_LAGRANGE) then 1860 1979 C … … 1906 2025 & s_parent,s_child,ds_parent,ds_child) 1907 2026 C 1908 elseif ( TYPEinterp .EQ. AGRIF_PPM ) then1909 Call ppm1D1910 & (parenttab,childtab,1911 & indmax-indmin+1,petab_child-pttab_child+1,1912 & s_parent,s_child,ds_parent,ds_child)1913 C1914 2027 endif 1915 2028 C -
trunk/AGRIF/AGRIF_FILES/modinterpbasic.F
r662 r779 37 37 C 38 38 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 40 47 41 48 CONTAINS … … 140 147 C 141 148 End Subroutine Linear1d 149 150 Subroutine Linear1dprecompute2d(np2, np,nc, 151 & s_parent,s_child,ds_parent,ds_child,dir) 152 C 153 CCC Description: 154 CCC Subroutine to compute 2D coefficient and index for a linear 1D interpolation 155 CCC on a child grid (vector y) from its parent grid (vector x). 156 C 157 CC Method: 158 C 159 C Declarations: 160 C 161 162 C 163 C Arguments 164 INTEGER :: np,nc,np2,dir 165 REAL :: s_parent,s_child,ds_parent,ds_child 166 C 167 C 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 174 C 175 C 176 177 coeffraf = nint(ds_parent/ds_child) 178 C 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 188 C 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 214 C 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 225 C 226 enddo 227 C 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 235 C 236 globind_parent_left = s_parent 237 & + (locind_parent_left - 1)*ds_parent 238 C 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 256 C 257 C 258 End Subroutine Linear1dprecompute2d 259 260 261 262 Subroutine Linear1dprecompute(np,nc, 263 & s_parent,s_child,ds_parent,ds_child) 264 C 265 CCC Description: 266 CCC Subroutine to compute 1D coefficient and index for a linear 1D interpolation 267 CCC on a child grid (vector y) from its parent grid (vector x). 268 C 269 CC Method: 270 C 271 C Declarations: 272 C 273 274 C 275 C Arguments 276 INTEGER :: np,nc 277 REAL :: s_parent,s_child,ds_parent,ds_child 278 C 279 C 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 284 C 285 C 286 287 coeffraf = nint(ds_parent/ds_child) 288 C 289 if (coeffraf == 1) then 290 C 291 return 292 C 293 endif 294 C 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 304 C 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 321 C 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 331 C 332 enddo 333 C 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 341 C 342 globind_parent_left = s_parent 343 & + (locind_parent_left - 1)*ds_parent 344 C 345 indparent(nc,1) = locind_parent_left 346 347 coeffparent(nc,1) = (globind_parent_left + ds_parent - ypos) 348 & * invds 349 endif 350 C 351 Return 352 C 353 C 354 End Subroutine Linear1dprecompute 355 356 Subroutine Linear1daftercompute(x,y,np,nc, 357 & s_parent,s_child,ds_parent,ds_child,dir) 358 C 359 CCC Description: 360 CCC Subroutine to do a linear 1D interpolation on a child grid (vector y) from 361 CCC its parent grid (vector x) using precomputed coefficient and index. 362 C 363 CC Method: 364 C 365 C Declarations: 366 C 367 368 C 369 C 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 374 C 375 C 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 380 C 381 C 382 383 !CDIR ALTCODE 384 !CDIR NODEP 385 do i = 1,nc 386 C 387 y(i)=coeffparent(i,dir)*x(indparent(i,dir))+ 388 & (1.-coeffparent(i,dir))*x(indparent(i,dir)+1) 389 C 390 enddo 391 C 392 Return 393 C 394 C 395 End Subroutine Linear1daftercompute 142 396 143 397 C … … 569 823 If (coeffraf == 1) Then 570 824 locind_parent_left = 1 + nint((s_child - s_parent)/ds_parent) 825 !CDIR ALTCODE 826 !CDIR SHORTLOOP 571 827 y(1:nc) = x(locind_parent_left:locind_parent_left+nc-1) 572 828 return … … 582 838 Allocate(tabtest4(-2*coeffraf:nc+2*coeffraf)) 583 839 ENDIF 584 ENDIF 840 ENDIF 841 585 842 ypos = s_child 586 843 C … … 595 852 C 596 853 854 !CDIR NOVECTOR 597 855 Do i=1,coeffraf 598 856 tabdiff2(i)=(real(i)-0.5)*invcoeffraf … … 602 860 tabdiff3(1) = (1./3.)*a 603 861 a=2.*a 862 !CDIR NOVECTOR 604 863 Do i=2,coeffraf 605 864 tabdiff3(i) = tabdiff3(i-1)+(real(i)-1)*a … … 621 880 C 622 881 C 882 !CDIR ALTCODE 883 !CDIR SHORTLOOP 623 884 Do i = nmin,nmax 624 885 slope(i) = x(i) - x(i-1) 625 886 Enddo 626 887 888 !CDIR ALTCODE 889 !CDIR SHORTLOOP 627 890 Do i = nmin+1,nmax-1 628 891 xl(i)= 0.5*(x(i-1)+x(i)) … … 631 894 C 632 895 C apply parabolic monotonicity 896 !CDIR ALTCODE 897 !CDIR SHORTLOOP 633 898 Do i = locind_parent_left,locind_parent_last 634 899 delta(i) = xl(i+1) - xl(i) … … 644 909 Do iparent = locind_parent_left,locind_parent_last 645 910 pos=1 911 !CDIR ALTCODE 912 !CDIR SHORTLOOP 646 913 Do jj = ipos - coeffraf/2+diffmod,ipos + coeffraf/2 647 914 C … … 653 920 ipos = ipos + coeffraf 654 921 C 655 End do 656 C 657 C 922 End do 923 C 924 C 925 !CDIR ALTCODE 926 !CDIR SHORTLOOP 658 927 y(1:nc)=tabtest4(1:nc) 659 928 660 929 Return 661 930 End Subroutine ppm1d 931 932 933 Subroutine ppm1dprecompute2d(np2,np,nc, 934 & s_parent,s_child,ds_parent,ds_child,dir) 935 C 936 CCC Description: 937 CCC Subroutine to compute 2D coefficients and index for a 1D interpolation 938 CCC using piecewise parabolic method 939 CC Method: 940 C 941 C Declarations: 942 C 943 Implicit none 944 C 945 C Arguments 946 Integer :: np2,np,nc,dir 947 C Real, Dimension(:),Allocatable :: ytemp 948 Real :: s_parent,s_child,ds_parent,ds_child 949 C 950 C 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 958 C 959 Real :: xrmin,xrmax,am3,s2,s1 960 Real, Dimension(np) :: xl,delta,a6,slope 961 INTEGER :: diffmod 962 REAL :: invcoeffraf 963 C 964 coeffraf = nint(ds_parent/ds_child) 965 C 966 invcoeffraf = ds_child/ds_parent 967 C 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 1007 C 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) 1012 C 1013 xpmin = s_parent + (locind_parent_left-1)*ds_parent 1014 i1 = 1+agrif_int((xpmin-s_child)/ds_child) 1015 C 1016 C 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 1043 C 1044 C 1045 diffmod = 0 1046 IF (mod(coeffraf,2) == 0) diffmod = 1 1047 C 1048 ipos = i1 1049 C 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 1059 C 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) 1077 C 1078 CCC Description: 1079 CCC Subroutine to compute coefficient and index for a 1D interpolation 1080 CCC using piecewise parabolic method 1081 CC Method: 1082 C 1083 C Declarations: 1084 C 1085 Implicit none 1086 C 1087 C Arguments 1088 Integer :: np,nc 1089 C Real, Dimension(:),Allocatable :: ytemp 1090 Real :: s_parent,s_child,ds_parent,ds_child 1091 C 1092 C 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 1098 C 1099 Real :: xrmin,xrmax,am3,s2,s1 1100 Real, Dimension(np) :: xl,delta,a6,slope 1101 C Real, Dimension(:),Allocatable :: diff,diff2,diff3 1102 INTEGER :: diffmod 1103 REAL :: invcoeffraf 1104 C 1105 coeffraf = nint(ds_parent/ds_child) 1106 C 1107 If (coeffraf == 1) Then 1108 return 1109 End If 1110 invcoeffraf = ds_child/ds_parent 1111 C 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 1125 C 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) 1130 C 1131 xpmin = s_parent + (locind_parent_left-1)*ds_parent 1132 i1 = 1+agrif_int((xpmin-s_child)/ds_child) 1133 C 1134 C 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 1160 C 1161 C 1162 diffmod = 0 1163 IF (mod(coeffraf,2) == 0) diffmod = 1 1164 C 1165 ipos = i1 1166 C 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 1177 C 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) 1185 C 1186 CCC Description: 1187 CCC Subroutine to do a 1D interpolation and apply monotonicity constraints 1188 CCC using piecewise parabolic method 1189 CCC on a child grid (vector y) from its parent grid (vector x). 1190 CC Method: Use precomputed coefficient and index 1191 C 1192 C Declarations: 1193 C 1194 Implicit none 1195 C 1196 C Arguments 1197 Integer :: np,nc,dir 1198 Real, INTENT(IN),Dimension(np) :: x 1199 Real, INTENT(OUT),Dimension(nc) :: y 1200 C Real, Dimension(:),Allocatable :: ytemp 1201 Real :: s_parent,s_child,ds_parent,ds_child 1202 C 1203 C 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 1209 C 1210 Real :: xrmin,xrmax,am3,s2,s1 1211 Real, Dimension(np) :: xl,delta,a6,slope 1212 INTEGER :: diffmod 1213 C 1214 C 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 662 1225 663 1226 C ************************************************************************** -
trunk/AGRIF/AGRIF_FILES/modmask.F
r662 r779 35 35 C 36 36 IMPLICIT NONE 37 Integer, Parameter :: MaxSearch = 537 Integer, Parameter :: MaxSearch = 3 38 38 C 39 39 CONTAINS … … 103 103 IF (tempP%var%array3(i0,j0,k0) 104 104 & == 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 112 c Call CalculNewValTempP((/i0,j0,k0/), 113 c & tempP,parent, 114 c & ppbtab,ppetab, 115 c & noraftab,nbdim) 116 109 117 ENDIF 110 118 enddo … … 118 126 IF (tempP%var%array4(i0,j0,k0,l0) 119 127 & == 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) 124 133 ENDIF 125 134 enddo … … 221 230 Valmax = min(Valmax,MaxSearch) 222 231 C 232 !CDIR NOVECTOR 223 233 imin = indic 234 !CDIR NOVECTOR 224 235 imax = indic 225 236 C … … 238 249 if (indic(iii).GT.ppbtab(iii)) then 239 250 251 !CDIR NOVECTOR 240 252 idecal = indic 241 253 idecal(iii) = idecal(iii)-1 … … 286 298 SELECT CASE(nbdim) 287 299 CASE (1) 300 !CDIR ALTCODE 301 !CDIR SHORTLOOP 288 302 do ii = imin(1),imax(1) 289 303 ValParent = parent%var%array1(ii) … … 296 310 CASE (2) 297 311 do jj = imin(2),imax(2) 312 !CDIR ALTCODE 313 !CDIR SHORTLOOP 298 314 do ii = imin(1),imax(1) 299 315 ValParent = parent%var%array2(ii,jj) … … 308 324 do kk = imin(3),imax(3) 309 325 do jj = imin(2),imax(2) 326 !CDIR ALTCODE 327 !CDIR SHORTLOOP 310 328 do ii = imin(1),imax(1) 311 329 ValParent = parent%var%array3(ii,jj,kk) … … 322 340 do kk = imin(3),imax(3) 323 341 do jj = imin(2),imax(2) 342 !CDIR ALTCODE 343 !CDIR SHORTLOOP 324 344 do ii = imin(1),imax(1) 325 345 ValParent = parent%var%array4(ii,jj,kk,ll) … … 338 358 do kk = imin(3),imax(3) 339 359 do jj = imin(2),imax(2) 360 !CDIR ALTCODE 361 !CDIR SHORTLOOP 340 362 do ii = imin(1),imax(1) 341 363 ValParent = parent%var%array5(ii,jj,kk,ll,mm) … … 356 378 do kk = imin(3),imax(3) 357 379 do jj = imin(2),imax(2) 380 !CDIR ALTCODE 381 !CDIR SHORTLOOP 358 382 do ii = imin(1),imax(1) 359 383 ValParent = parent%var%array6(ii,jj,kk,ll,mm,nn) … … 412 436 C 413 437 C 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) 443 C 444 CCC Description: 445 CCC Subroutine called in the procedure Agrif_InterpnD to recalculate the value 446 CCC of the parent grid variable when this one is equal to Agrif_SpecialValue. 447 C 448 C Declarations: 449 C 450 451 C 452 C 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 458 C 459 C 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 463 C 464 C 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 473 C 474 C Local arrays 475 C 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 484 C 485 Valmax = min(Valmax,MaxSearch) 486 C 487 !CDIR NOVECTOR 488 imin = indic 489 !CDIR NOVECTOR 490 imax = indic 491 492 !CDIR NOVECTOR 493 idecal = indic 494 C 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 516 C 517 Existunmasked = .FALSE. 518 C 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 530 C 531 C 532 if (.Not.Existunmasked) return 533 C 534 i = 1 535 C 536 do While(i <= ValMax) 537 C 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 545 C 546 Res = 0. 547 Nbvals = 0 548 C 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 561 C 562 C 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 571 C 572 End Subroutine CalculNewValTempP3D 573 574 Subroutine CalculNewValTempP4D(indic, 575 & tempP,parent,ppbtab, 576 & ppetab,noraftab,MaxSearch,Agrif_SpecialValue) 577 C 578 CCC Description: 579 CCC Subroutine called in the procedure Agrif_InterpnD to recalculate the value 580 CCC of the parent grid variable when this one is equal to Agrif_SpecialValue. 581 C 582 C Declarations: 583 C 584 585 C 586 C 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 593 C 594 C 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 599 C 600 C 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 608 C 609 C Local arrays 610 C 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 618 C 619 Valmax = min(Valmax,MaxSearch) 620 C 621 imin = indic 622 imax = indic 623 C 624 i = 1 625 firsttest = .TRUE. 626 idecal = indic 627 628 C 629 do While(i <= ValMax) 630 C 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 653 C 654 Res = 0. 655 Nbvals = 0 656 C 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 670 C 671 C 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 688 C 689 End Subroutine CalculNewValTempP4D -
trunk/AGRIF/AGRIF_FILES/modmpp.F
r774 r779 146 146 sendtoproc(k) = .true. 147 147 C 148 !CDIR SHORTLOOP 148 149 do i = 1,nbdim 149 150 C … … 228 229 datasize = 1 229 230 C 231 !CDIR SHORTLOOP 230 232 do i = 1,nbdim 231 233 C … … 305 307 datasize = 1 306 308 C 309 !CDIR SHORTLOOP 307 310 do i = 1,nbdim 308 311 C -
trunk/AGRIF/AGRIF_FILES/modupdate.F
r662 r779 42 42 C 43 43 IMPLICIT NONE 44 logical, private :: precomputedone(7) = .FALSE. 44 45 C 45 46 CONTAINS … … 1301 1302 1302 1303 if ( nbdim .EQ. 1 ) then 1304 tempP%var%array1 = 0. 1303 1305 Call Agrif_Update_1D_recursive(TypeUpdate, 1304 1306 & tempP%var%array1,tempCextend%var%array1, … … 1822 1824 REAL, DIMENSION(indmin(1):indmax(1), 1823 1825 & 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 1824 1830 INTEGER :: i,j 1825 1831 INTEGER :: coeffraf,locind_child_left 1826 1832 C 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 1827 1866 do j = pttab_child(nbdim),petab_child(nbdim) 1828 1867 C 1868 !---CDIR NEXPAND 1829 1869 Call Agrif_Update_1D_recursive(TypeUpdate(1:nbdim-1), 1830 1870 & tabtemp(:,j), … … 1836 1876 C 1837 1877 enddo 1838 1878 ENDIF 1879 tabtemp_trsp = TRANSPOSE(tabtemp) 1880 coeffraf = nint(ds_parent(nbdim)/ds_child(nbdim)) 1881 1882 !---CDIR NEXPAND 1839 1883 Call Agrif_Compute_nbdim_update(s_parent(nbdim),s_child(nbdim), 1840 1884 & ds_parent(nbdim),ds_child(nbdim),coeffraf,locind_child_left) 1841 1885 C 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 1842 1919 do i = indmin(1),indmax(1) 1843 1920 C 1921 !---CDIR NEXPAND 1844 1922 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), 1847 1925 & indmin(nbdim),indmax(nbdim), 1848 1926 & pttab_child(nbdim),petab_child(nbdim), … … 1852 1930 C 1853 1931 enddo 1932 1933 ENDIF 1934 1935 tempP = TRANSPOSE(tempP_trsp) 1854 1936 C 1855 1937 Return … … 1902 1984 C 1903 1985 C 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 1904 2029 do k = pttab_child(nbdim),petab_child(nbdim) 1905 2030 C … … 1914 2039 enddo 1915 2040 C 2041 precomputedone(1) = .FALSE. 2042 precomputedone(2) = .FALSE. 2043 coeffraf = nint ( ds_parent(3) / ds_child(3) ) 2044 1916 2045 Call Agrif_Compute_nbdim_update(s_parent(nbdim),s_child(nbdim), 1917 2046 & ds_parent(nbdim),ds_child(nbdim),coeffraf,locind_child_left) … … 1930 2059 1931 2060 ELSE 2061 tempP = 0. 1932 2062 C 1933 2063 do j = indmin(2),indmax(2) … … 2021 2151 & ds_parent(nbdim),ds_child(nbdim),coeffraf,locind_child_left) 2022 2152 C 2153 tempP = 0. 2154 2023 2155 do k = indmin(3),indmax(3) 2024 2156 C … … 2120 2252 Call Agrif_Compute_nbdim_update(s_parent(nbdim),s_child(nbdim), 2121 2253 & ds_parent(nbdim),ds_child(nbdim),coeffraf,locind_child_left) 2254 tempP = 0. 2122 2255 C 2123 2256 do l = indmin(4),indmax(4) … … 2232 2365 & ds_parent(nbdim),ds_child(nbdim),coeffraf,locind_child_left) 2233 2366 C 2367 tempP = 0. 2368 2234 2369 do m = indmin(5),indmax(5) 2235 2370 do l = indmin(4),indmax(4) … … 2307 2442 C 2308 2443 elseif (TypeUpdate == AGRIF_Update_average) then 2309 C 2444 C 2310 2445 Call average1D 2311 2446 & (parenttab,childtab, 2312 2447 & 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) 2314 2449 C 2315 2450 elseif (TypeUpdate == AGRIF_Update_full_weighting) then … … 2421 2556 #endif 2422 2557 2423 End Module Agrif_Update2424 2558 2425 2559 C ************************************************************************** … … 2443 2577 C 2444 2578 C Arguments 2445 USE Agrif_Update2446 2579 INTEGER :: nbdim 2447 2580 INTEGER, DIMENSION(nbdim) :: TypeUpdate ! TYPE of update (copy or average) … … 2466 2599 & ds_parent(nbdim),ds_child(nbdim), 2467 2600 & coeffraf,locind_child_left) 2601 2468 2602 C 2469 2603 Return … … 2471 2605 C 2472 2606 End Subroutine Agrif_Update_1D_recursive 2607 2608 End Module Agrif_Update -
trunk/AGRIF/AGRIF_FILES/modupdatebasic.F
r662 r779 37 37 38 38 IMPLICIT NONE 39 40 integer,dimension(:,:),allocatable :: indchildcopy 41 integer,dimension(:,:),allocatable :: indchildaverage 39 42 C 40 43 … … 78 81 locind_child_left = 1 + nint((s_parent - s_child)/ds_child) 79 82 C 83 !CDIR ALTCODE 80 84 x(1:np) = y(locind_child_left:locind_child_left+np-1) 81 85 C … … 87 91 locind_child_left = 1 + nint((s_parent - s_child)/ds_child) 88 92 93 !CDIR ALTCODE 89 94 do i = 1,np 90 95 C … … 100 105 C 101 106 End Subroutine agrif_copy1d 102 C 103 C 107 108 C ************************************************************************** 109 CCC Subroutine Copy1dprecompute 110 C ************************************************************************** 111 C 112 Subroutine copy1dprecompute2d(nc2,np,nc, 113 & s_parent,s_child,ds_parent,ds_child,dir) 114 C 115 CCC Description: 116 CCC Subroutine to precompute index for a copy on a parent grid (vector x) from 117 CCC its child grid (vector y). 118 C 119 CC Method: 120 C 121 C Declarations: 122 C 123 124 C 125 C Arguments 126 INTEGER :: nc2,np,nc 127 INTEGER :: dir 128 REAL :: s_parent,s_child 129 REAL :: ds_parent,ds_child 130 C 131 C Local variables 132 INTEGER,DIMENSION(:,:),ALLOCATABLE :: indchildcopy_tmp 133 INTEGER :: i,locind_child_left,coeffraf 134 C 135 C 136 coeffraf = nint(ds_parent/ds_child) 137 C 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 157 C 158 indchildcopy(i,dir) = locind_child_left 159 locind_child_left = locind_child_left + coeffraf 160 C 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 168 C 169 Return 170 C 171 C 172 End Subroutine copy1dprecompute2d 173 174 C 175 C 176 C ************************************************************************** 177 CCC Subroutine Copy1dprecompute 178 C ************************************************************************** 179 C 180 Subroutine copy1dprecompute(np,nc, 181 & s_parent,s_child,ds_parent,ds_child,dir) 182 C 183 CCC Description: 184 CCC Subroutine to precompute index for a copy on a parent grid (vector x) from 185 CCC its child grid (vector y). 186 C 187 CC Method: 188 C 189 C Declarations: 190 C 191 192 C 193 C Arguments 194 INTEGER :: np,nc 195 INTEGER :: dir 196 REAL :: s_parent,s_child 197 REAL :: ds_parent,ds_child 198 C 199 C Local variables 200 INTEGER :: i,locind_child_left,coeffraf 201 C 202 C 203 coeffraf = nint(ds_parent/ds_child) 204 C 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 218 C 219 indchildcopy(i,1) = locind_child_left 220 locind_child_left = locind_child_left + coeffraf 221 C 222 enddo 223 224 C 225 Return 226 C 227 C 228 End Subroutine copy1dprecompute 229 C 230 C 231 C ************************************************************************** 232 CCC Subroutine Copy1daftercompute 233 C ************************************************************************** 234 C 235 Subroutine copy1daftercompute(x,y,np,nc, 236 & s_parent,s_child,ds_parent,ds_child,dir) 237 C 238 CCC Description: 239 CCC Subroutine to do a copy on a parent grid (vector x) from its child grid 240 CCC (vector y) using precomputed index. 241 C 242 CC Method: 243 C 244 C Declarations: 245 C 246 247 C 248 C 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 255 C 256 C Local variables 257 INTEGER :: i 258 C 259 C 260 261 !CDIR ALTCODE 262 do i = 1,np 263 C 264 x(i) = y(indchildcopy(i,dir)) 265 C 266 enddo 267 268 C 269 Return 270 C 271 C 272 End Subroutine copy1daftercompute 273 274 104 275 C 105 276 C ************************************************************************** … … 123 294 C Local variables 124 295 INTEGER :: i,locind_child_left,coeffraf,ii 125 REAL :: xpos 296 REAL :: xpos, invcoeffraf 126 297 INTEGER :: nbnonnuls 127 298 INTEGER :: diffmod … … 129 300 C 130 301 coeffraf = nint(ds_parent/ds_child) 302 invcoeffraf = 1./coeffraf 131 303 C 132 304 if (coeffraf == 1) then … … 149 321 150 322 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 323 C 324 IF (Agrif_UseSpecialValueInUpdate) THEN 325 326 do i = 1,np 162 327 nbnonnuls = 0 328 !CDIR NOVECTOR 163 329 Do ii = -coeffraf/2+locind_child_left+diffmod, 164 330 & coeffraf/2+locind_child_left 165 331 C 166 IF (Agrif_UseSpecialValueInUpdate) THEN 167 IF (y(ii) .NE. Agrif_SpecialValueFineGrid) THEN 332 IF (y(ii) .NE. Agrif_SpecialValueFineGrid) THEN 168 333 nbnonnuls = nbnonnuls + 1 169 334 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 352 C 172 353 x(i) = x(i) + y(ii) 173 ENDIF174 354 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 188 356 locind_child_left = locind_child_left + coeffraf 189 C 190 enddo 191 C 357 enddo 358 359 ENDIF 360 192 361 Return 193 362 C 194 363 C 195 364 End Subroutine average1d 365 366 Subroutine average1dprecompute2d(nc2,np,nc, 367 & s_parent,s_child,ds_parent,ds_child,dir) 368 C 369 CCC Description: 370 CCC Subroutine to do an update by average on a parent grid (vector x)from its 371 CCC child grid (vector y). 372 C 373 C Arguments 374 INTEGER :: nc2,np,nc,dir 375 REAL :: s_parent,s_child 376 REAL :: ds_parent,ds_child 377 C 378 C 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 384 C 385 C 386 coeffraf = nint(ds_parent/ds_child) 387 C 388 xpos = s_parent 389 C 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 422 C 423 C 424 End Subroutine average1dprecompute2d 425 426 427 Subroutine average1dprecompute(np,nc, 428 & s_parent,s_child,ds_parent,ds_child) 429 C 430 CCC Description: 431 CCC Subroutine to do an update by average on a parent grid (vector x)from its 432 CCC child grid (vector y). 433 C 434 C Arguments 435 INTEGER :: np,nc 436 REAL :: s_parent,s_child 437 REAL :: ds_parent,ds_child 438 C 439 C Local variables 440 INTEGER :: i,locind_child_left,coeffraf,ii 441 REAL :: xpos, invcoeffraf 442 INTEGER :: nbnonnuls 443 INTEGER :: diffmod 444 C 445 C 446 coeffraf = nint(ds_parent/ds_child) 447 448 C 449 if (coeffraf == 1) then 450 C 451 return 452 C 453 endif 454 455 C 456 xpos = s_parent 457 C 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 475 C 476 indchildaverage(i,1)=-coeffraf/2+locind_child_left+diffmod 477 locind_child_left = locind_child_left + coeffraf 478 enddo 479 480 Return 481 C 482 C 483 End Subroutine average1dprecompute 484 485 Subroutine average1daftercompute(x,y,np,nc, 486 & s_parent,s_child,ds_parent,ds_child,dir) 487 C 488 CCC Description: 489 CCC Subroutine to do an update by average on a parent grid (vector x)from its 490 CCC child grid (vector y). 491 C 492 C 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 498 C 499 C 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/) 506 C 507 C 508 coeffraf = nint(ds_parent/ds_child) 509 invcoeffraf = 1./coeffraf 510 C 511 C 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 540 C 541 C 542 End Subroutine average1daftercompute 196 543 C 197 544 C
Note: See TracChangeset
for help on using the changeset viewer.