|
Oasis3 4.0.2
|
|
Oasis3 4.0.2
|
00001 C**** 00002 C ***************************** 00003 C * OASIS MODULE - LEVEL ? * 00004 C * ------------- ------- * 00005 C ***************************** 00006 C 00007 C**** remap_bilinear - calculate bilinear remapping 00008 C 00009 C Purpose: 00010 C ------- 00011 C Adaptation of SCRIP 1.4 remap_bicubic module to calculate 00012 C bicubic remapping. 00013 C 00014 C** Interface: 00015 C --------- 00016 C *CALL* *remap_bicub** 00017 C 00018 C Input: 00019 C ----- 00020 C 00021 C Output: 00022 C ------ 00023 C 00024 C History: 00025 C ------- 00026 C Version Programmer Date Description 00027 C ------- ---------- ---- ----------- 00028 C 2.5 S. Valcke 2002/09 Treatment for masked point 00029 C 00030 C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 00031 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00032 ! 00033 ! this module contains necessary routines for performing an 00034 ! bicubic interpolation. 00035 ! 00036 !----------------------------------------------------------------------- 00037 ! 00038 ! CVS:$Id: remap_bicubic.f 3500 2012-03-22 14:58:30Z coquart $ 00039 ! 00040 ! Copyright (c) 1997, 1998 the Regents of the University of 00041 ! California. 00042 ! 00043 ! This software and ancillary information (herein called software) 00044 ! called SCRIP is made available under the terms described here. 00045 ! The software has been approved for release with associated 00046 ! LA-CC Number 98-45. 00047 ! 00048 ! Unless otherwise indicated, this software has been authored 00049 ! by an employee or employees of the University of California, 00050 ! operator of the Los Alamos National Laboratory under Contract 00051 ! No. W-7405-ENG-36 with the U.S. Department of Energy. The U.S. 00052 ! Government has rights to use, reproduce, and distribute this 00053 ! software. The public may copy and use this software without 00054 ! charge, provided that this Notice and any statement of authorship 00055 ! are reproduced on all copies. Neither the Government nor the 00056 ! University makes any warranty, express or implied, or assumes 00057 ! any liability or responsibility for the use of this software. 00058 ! 00059 ! If software is modified to produce derivative works, such modified 00060 ! software should be clearly marked, so as not to confuse it with 00061 ! the version available from Los Alamos National Laboratory. 00062 ! 00063 !*********************************************************************** 00064 00065 module remap_bicubic 00066 00067 !----------------------------------------------------------------------- 00068 00069 use kinds_mod ! defines common data types 00070 use constants ! defines common constants 00071 use grids ! module containing grid info 00072 use remap_vars ! module containing remap info 00073 00074 implicit NONE 00075 00076 !----------------------------------------------------------------------- 00077 00078 integer (kind=int_kind), parameter :: 00079 & max_iter = 100 ! max iteration count for i,j iteration 00080 00081 real (kind=dbl_kind), parameter :: 00082 & converge = epsilon(1.0_dbl_kind) ! convergence criterion 00083 00084 !*********************************************************************** 00085 00086 contains 00087 00088 !*********************************************************************** 00089 00090 subroutine remap_bicub (lextrapdone) 00091 00092 !----------------------------------------------------------------------- 00093 ! 00094 ! this routine computes the weights for a bicubic interpolation. 00095 ! 00096 !----------------------------------------------------------------------- 00097 00098 LOGICAL :: lextrapdone, ! logical, true if EXTRAP done on field 00099 & ll_nnei ! true if extra search is done 00100 !----------------------------------------------------------------------- 00101 ! 00102 ! local variables 00103 ! 00104 !----------------------------------------------------------------------- 00105 00106 integer (kind=int_kind) :: n,icount,i, 00107 & dst_add, ! destination address 00108 & iter, ! iteration counter 00109 & nmap ! index of current map being computed 00110 00111 integer (kind=int_kind), dimension(4) :: 00112 & src_add ! address for the four source points 00113 00114 real (kind=dbl_kind), dimension(4) :: 00115 & src_lats, ! latitudes of four corners 00116 & src_lons ! longitudes of four corners 00117 00118 real (kind=dbl_kind), dimension(4,4) :: 00119 & wgts ! bicubic weights for four corners 00120 00121 real (kind=dbl_kind) :: 00122 & plat, plon, ! lat/lon coords of destination point 00123 & iguess, jguess, ! current guess for bicubic coordinate 00124 & deli, delj, ! corrections to i,j 00125 & dth1, dth2, dth3, ! some latitude differences 00126 & dph1, dph2, dph3, ! some longitude differences 00127 & dthp, dphp, ! difference between point and sw corner 00128 & mat1, mat2, mat3, mat4, ! matrix elements 00129 & determinant, ! matrix determinant 00130 & sum_wgts ! sum of weights for normalization 00131 00132 real (kind=dbl_kind) :: 00133 & coslat_dst, sinlat_dst, coslon_dst, sinlon_dst, 00134 & dist_min, distance, ! for computing dist-weighted avg 00135 & src_latsnn, arg 00136 00137 integer (kind=int_kind) :: min_add, max_add, srch_add, src_addnn 00138 00139 !----------------------------------------------------------------------- 00140 ! 00141 ! compute mappings from grid1 to grid2 00142 ! 00143 !----------------------------------------------------------------------- 00144 ! 00145 IF (nlogprt .GE. 2) THEN 00146 WRITE (UNIT = nulou,FMT = *)' ' 00147 WRITE (UNIT = nulou,FMT = *)'Entering routine remap_bicub' 00148 CALL FLUSH(nulou) 00149 ENDIF 00150 ! 00151 ll_nnei = .true. 00152 nmap = 1 00153 if (grid1_rank /= 2) then 00154 stop 'Can not do bicubic interpolation when grid_rank /= 2' 00155 endif 00156 00157 !*** 00158 !*** loop over destination grid 00159 !*** 00160 00161 grid_loop1: do dst_add = 1, grid2_size 00162 00163 if (.not. grid2_mask(dst_add)) cycle grid_loop1 00164 00165 plat = grid2_center_lat(dst_add) 00166 plon = grid2_center_lon(dst_add) 00167 00168 !*** 00169 !*** find nearest square of grid points on source grid 00170 !*** 00171 00172 call grid_search_bicub(src_add, src_lats, src_lons, 00173 & min_add, max_add, 00174 & plat, plon, grid1_dims, 00175 & grid1_center_lat, grid1_center_lon, 00176 & grid1_bound_box, bin_addr1, 00177 & lextrapdone) 00178 00179 if (src_add(1) > 0) then 00180 00181 !*** 00182 !*** if the 4 surrounding points have been found and are 00183 !*** non-masked, do bicubic interpolation 00184 !*** 00185 00186 grid2_frac(dst_add) = one 00187 00188 dth1 = src_lats(2) - src_lats(1) 00189 dth2 = src_lats(4) - src_lats(1) 00190 dth3 = src_lats(3) - src_lats(2) - dth2 00191 00192 dph1 = src_lons(2) - src_lons(1) 00193 dph2 = src_lons(4) - src_lons(1) 00194 dph3 = src_lons(3) - src_lons(2) 00195 00196 if (dph1 > three*pih) dph1 = dph1 - pi2 00197 if (dph2 > three*pih) dph2 = dph2 - pi2 00198 if (dph3 > three*pih) dph3 = dph3 - pi2 00199 if (dph1 < -three*pih) dph1 = dph1 + pi2 00200 if (dph2 < -three*pih) dph2 = dph2 + pi2 00201 if (dph3 < -three*pih) dph3 = dph3 + pi2 00202 00203 dph3 = dph3 - dph2 00204 00205 iguess = half 00206 jguess = half 00207 00208 iter_loop1: do iter=1,max_iter 00209 00210 dthp = plat - src_lats(1) - dth1*iguess - 00211 & dth2*jguess - dth3*iguess*jguess 00212 dphp = plon - src_lons(1) 00213 00214 if (dphp > three*pih) dphp = dphp - pi2 00215 if (dphp < -three*pih) dphp = dphp + pi2 00216 00217 dphp = dphp - dph1*iguess - dph2*jguess - 00218 & dph3*iguess*jguess 00219 00220 mat1 = dth1 + dth3*jguess 00221 mat2 = dth2 + dth3*iguess 00222 mat3 = dph1 + dph3*jguess 00223 mat4 = dph2 + dph3*iguess 00224 00225 determinant = mat1*mat4 - mat2*mat3 00226 00227 deli = (dthp*mat4 - mat2*dphp)/determinant 00228 delj = (mat1*dphp - dthp*mat3)/determinant 00229 00230 if (abs(deli) < converge .and. 00231 & abs(delj) < converge) exit iter_loop1 00232 00233 iguess = iguess + deli 00234 jguess = jguess + delj 00235 00236 end do iter_loop1 00237 00238 if (iter <= max_iter) then 00239 00240 !*** 00241 !*** successfully found i,j - compute weights 00242 !*** 00243 00244 wgts(1,1) = (one - jguess**2*(three-two*jguess))* 00245 & (one - iguess**2*(three-two*iguess)) 00246 wgts(1,2) = (one - jguess**2*(three-two*jguess))* 00247 & iguess**2*(three-two*iguess) 00248 wgts(1,3) = jguess**2*(three-two*jguess)* 00249 & iguess**2*(three-two*iguess) 00250 wgts(1,4) = jguess**2*(three-two*jguess)* 00251 & (one - iguess**2*(three-two*iguess)) 00252 wgts(2,1) = (one - jguess**2*(three-two*jguess))* 00253 & iguess*(iguess-one)**2 00254 wgts(2,2) = (one - jguess**2*(three-two*jguess))* 00255 & iguess**2*(iguess-one) 00256 wgts(2,3) = jguess**2*(three-two*jguess)* 00257 & iguess**2*(iguess-one) 00258 wgts(2,4) = jguess**2*(three-two*jguess)* 00259 & iguess*(iguess-one)**2 00260 wgts(3,1) = jguess*(jguess-one)**2* 00261 & (one - iguess**2*(three-two*iguess)) 00262 wgts(3,2) = jguess*(jguess-one)**2* 00263 & iguess**2*(three-two*iguess) 00264 wgts(3,3) = jguess**2*(jguess-one)* 00265 & iguess**2*(three-two*iguess) 00266 wgts(3,4) = jguess**2*(jguess-one)* 00267 & (one - iguess**2*(three-two*iguess)) 00268 wgts(4,1) = iguess*(iguess-one)**2* 00269 & jguess*(jguess-one)**2 00270 wgts(4,2) = iguess**2*(iguess-one)* 00271 & jguess*(jguess-one)**2 00272 wgts(4,3) = iguess**2*(iguess-one)* 00273 & jguess**2*(jguess-one) 00274 wgts(4,4) = iguess*(iguess-one)**2* 00275 & jguess**2*(jguess-one) 00276 00277 call store_link_bicub(dst_add, src_add, wgts, nmap) 00278 00279 ELSE 00280 WRITE (UNIT = nulou,FMT = *) 00281 & 'Iteration for i,j exceed max iteration count' 00282 STOP 00283 endif 00284 00285 else if (src_add(1) < 0) THEN 00286 00287 !*** 00288 !*** Search for bicubic failed or at least one of the 4 00289 !*** neighbours was masked. Do distance-weighted average using 00290 !*** the non-masked points among the 4 closest ones. 00291 !*** 00292 00293 00294 IF (nlogprt .eq. 2) then 00295 WRITE(nulou,*) ' ' 00296 WRITE(nulou,*) 00297 & 'WARNING: Cannot make bicubic interpolation for target point' 00298 & ,dst_add 00299 WRITE(nulou,*) 00300 & 'Using non-masked points among 4 nearest neighbors.' 00301 WRITE(nulou,*) ' ' 00302 ENDIF 00303 00304 !*** 00305 !*** Find the 4 closest points 00306 !*** 00307 coslat_dst = cos(plat) 00308 sinlat_dst = sin(plat) 00309 coslon_dst = cos(plon) 00310 sinlon_dst = sin(plon) 00311 src_add = 0 00312 dist_min = bignum 00313 src_lats = bignum 00314 do srch_add = min_add,max_add 00315 arg = coslat_dst*cos(grid1_center_lat(srch_add))* 00316 & (coslon_dst*cos(grid1_center_lon(srch_add)) + 00317 & sinlon_dst*sin(grid1_center_lon(srch_add)))+ 00318 & sinlat_dst*sin(grid1_center_lat(srch_add)) 00319 IF (arg < -1.0d0) THEN 00320 arg = -1.0d0 00321 ELSE IF (arg > 1.0d0) THEN 00322 arg = 1.0d0 00323 END IF 00324 distance=acos(arg) 00325 if (distance < dist_min) then 00326 sort_loop: do n=1,4 00327 if (distance < src_lats(n)) then 00328 do i=4,n+1,-1 00329 src_add (i) = src_add (i-1) 00330 src_lats(i) = src_lats(i-1) 00331 end do 00332 src_add (n) = srch_add 00333 src_lats(n) = distance 00334 dist_min = src_lats(4) 00335 exit sort_loop 00336 endif 00337 end do sort_loop 00338 endif 00339 end do 00340 00341 src_lons = one/(src_lats + tiny) 00342 distance = sum(src_lons) 00343 src_lats = src_lons/distance 00344 00345 !*** 00346 !*** Among 4 closest points, keep only the non-masked ones 00347 !*** 00348 00349 icount = 0 00350 do n=1,4 00351 if (grid1_mask(src_add(n)) .or. 00352 & (.not. grid1_mask(src_add(n)) .and. lextrapdone)) then 00353 icount = icount + 1 00354 else 00355 src_lats(n) = zero 00356 endif 00357 end do 00358 00359 if (icount > 0) then 00360 !*** renormalize weights 00361 sum_wgts = sum(src_lats) 00362 wgts(1,1) = src_lats(1)/sum_wgts 00363 wgts(1,2) = src_lats(2)/sum_wgts 00364 wgts(1,3) = src_lats(3)/sum_wgts 00365 wgts(1,4) = src_lats(4)/sum_wgts 00366 wgts(2,:) = 0. 00367 wgts(3,:) = 0. 00368 wgts(4,:) = 0. 00369 00370 grid2_frac(dst_add) = one 00371 call store_link_bicub(dst_add, src_add, wgts, nmap) 00372 ELSE 00373 IF (ll_nnei .eqv. .true. ) then 00374 IF (NLOGPRT .GE. 2) THEN 00375 WRITE(nulou,*) ' ' 00376 WRITE(nulou,*) 00377 ' Using the nearest non-masked neighbour & as all 4 surrounding source points are masked for target point', 00378 & dst_add 00379 WRITE(nulou,*) 'with longitude and latitude', 00380 & plon, plat 00381 ENDIF 00382 src_latsnn = bignum 00383 !cdir novector 00384 do srch_add = min_add,max_add 00385 if (grid1_mask(srch_add) .or. 00386 & (.not. grid1_mask(srch_add) .and. lextrapdone)) THEN 00387 arg = coslat_dst*cos(grid1_center_lat(srch_add))* 00388 & (coslon_dst*cos(grid1_center_lon(srch_add)) + 00389 & sinlon_dst*sin(grid1_center_lon(srch_add)))+ 00390 & sinlat_dst*sin(grid1_center_lat(srch_add)) 00391 IF (arg < -1.0d0) THEN 00392 arg = -1.0d0 00393 ELSE IF (arg > 1.0d0) THEN 00394 arg = 1.0d0 00395 END IF 00396 distance=acos(arg) 00397 00398 if (distance < src_latsnn) then 00399 src_addnn = srch_add 00400 src_latsnn = distance 00401 endif 00402 endif 00403 end DO 00404 IF (NLOGPRT .GE. 2) THEN 00405 WRITE(nulou,*) ' ' 00406 WRITE(nulou,*) 00407 & 'Nearest non masked neighbour is source point ', 00408 & src_addnn 00409 WRITE(nulou,*) 'with longitude and latitude', 00410 & grid1_center_lon(src_addnn), grid1_center_lat(src_addnn) 00411 ENDIF 00412 ! 00413 wgts(1,1) = 1. 00414 wgts(1,2:4) = 0. 00415 wgts(2,:) = 0. 00416 wgts(3,:) = 0. 00417 wgts(4,:) = 0. 00418 src_add(1) = src_addnn 00419 src_add(2) = 0 00420 src_add(3) = 0 00421 src_add(4) = 0 00422 00423 grid2_frac(dst_add) = one 00424 call store_link_bicub(dst_add, src_add, wgts, nmap) 00425 00426 ENDIF 00427 endif 00428 ENDIF 00429 end do grid_loop1 00430 ! 00431 IF (nlogprt .GE. 2) THEN 00432 WRITE (UNIT = nulou,FMT = *)' ' 00433 WRITE (UNIT = nulou,FMT = *)'Leaving routine remap_bicub' 00434 CALL FLUSH(nulou) 00435 ENDIF 00436 ! 00437 end subroutine remap_bicub 00438 00439 !*********************************************************************** 00440 00441 subroutine grid_search_bicub(src_add, src_lats, src_lons, 00442 & min_add, max_add, 00443 & plat, plon, src_grid_dims, 00444 & src_center_lat, src_center_lon, 00445 & src_bound_box, 00446 & src_bin_add, lextrapdone) 00447 00448 !----------------------------------------------------------------------- 00449 ! 00450 ! this routine finds the location of the search point plat, plon 00451 ! in the source grid and returns the corners needed for a bicubic 00452 ! interpolation. 00453 ! 00454 !----------------------------------------------------------------------- 00455 00456 !----------------------------------------------------------------------- 00457 ! 00458 ! output variables 00459 ! 00460 !----------------------------------------------------------------------- 00461 00462 integer (kind=int_kind), dimension(4), intent(out) :: 00463 & src_add ! address of each corner point enclosing P 00464 00465 real (kind=dbl_kind), dimension(4), intent(out) :: 00466 & src_lats, ! latitudes of the four corner points 00467 & src_lons ! longitudes of the four corner points 00468 00469 integer (kind=int_kind) :: min_add, max_add 00470 00471 !----------------------------------------------------------------------- 00472 ! 00473 ! input variables 00474 ! 00475 !----------------------------------------------------------------------- 00476 00477 real (kind=dbl_kind), intent(in) :: 00478 & plat, ! latitude of the search point 00479 & plon ! longitude of the search point 00480 00481 integer (kind=int_kind), dimension(2), intent(in) :: 00482 & src_grid_dims ! size of each src grid dimension 00483 00484 real (kind=dbl_kind), dimension(:), intent(in) :: 00485 & src_center_lat, ! latitude of each src grid center 00486 & src_center_lon ! longitude of each src grid center 00487 00488 real (kind=dbl_kind), dimension(:,:), intent(in) :: 00489 & src_bound_box ! bounding box for src grid search 00490 00491 integer (kind=int_kind), dimension(:,:), intent(in) :: 00492 & src_bin_add ! search bins for restricting 00493 00494 LOGICAL :: 00495 & lextrapdone ! logical, true if EXTRAP done on field 00496 00497 !----------------------------------------------------------------------- 00498 ! 00499 ! local variables 00500 ! 00501 !----------------------------------------------------------------------- 00502 00503 integer (kind=int_kind) :: n, next_n, srch_add, ni, ! dummy indices 00504 & nx, ny, ntotmask, ! dimensions of src grid 00505 & i, j, jp1, ip1, n_add, e_add, ne_add ! addresses 00506 00507 real (kind=dbl_kind) :: ! vectors for cross-product check 00508 & vec1_lat, vec1_lon, 00509 & vec2_lat, vec2_lon, cross_product, cross_product_last 00510 00511 !----------------------------------------------------------------------- 00512 ! 00513 ! restrict search first using search bins. 00514 ! 00515 !----------------------------------------------------------------------- 00516 00517 src_add = 0 00518 00519 min_add = size(src_center_lat) 00520 max_add = 1 00521 do n=1,num_srch_bins 00522 if (plat >= bin_lats(1,n) .and. plat <= bin_lats(2,n) .and. 00523 & plon >= bin_lons(1,n) .and. plon <= bin_lons(2,n)) then 00524 min_add = min(min_add, src_bin_add(1,n)) 00525 max_add = max(max_add, src_bin_add(2,n)) 00526 endif 00527 end do 00528 00529 !----------------------------------------------------------------------- 00530 ! 00531 ! now perform a more detailed search 00532 ! 00533 !----------------------------------------------------------------------- 00534 00535 nx = src_grid_dims(1) 00536 ny = src_grid_dims(2) 00537 00538 srch_loop: do srch_add = min_add,max_add 00539 00540 !*** first check bounding box 00541 00542 if (plat <= src_bound_box(2,srch_add) .and. 00543 & plat >= src_bound_box(1,srch_add) .and. 00544 & plon <= src_bound_box(4,srch_add) .and. 00545 & plon >= src_bound_box(3,srch_add)) then 00546 00547 !*** 00548 !*** we are within bounding box so get really serious 00549 !*** 00550 00551 !*** find N,S and NE points to this grid point 00552 00553 j = (srch_add - 1)/nx +1 00554 i = srch_add - (j-1)*nx 00555 00556 !*** find ip1 00557 !*** Note: I do not want to take into account the number 00558 !*** of overlapping grid points, as the underlying cell 00559 !*** will be found in all cases if the grid overlaps. 00560 00561 if (i < nx) then 00562 ip1 = i + 1 00563 else 00564 ip1 = 1 00565 endif 00566 00567 if (j < ny) then 00568 jp1 = j+1 00569 else 00570 jp1 = 1 00571 endif 00572 00573 n_add = (jp1 - 1)*nx + i 00574 e_add = (j - 1)*nx + ip1 00575 ne_add = (jp1 - 1)*nx + ip1 00576 00577 src_lats(1) = src_center_lat(srch_add) 00578 src_lats(2) = src_center_lat(e_add) 00579 src_lats(3) = src_center_lat(ne_add) 00580 src_lats(4) = src_center_lat(n_add) 00581 00582 src_lons(1) = src_center_lon(srch_add) 00583 src_lons(2) = src_center_lon(e_add) 00584 src_lons(3) = src_center_lon(ne_add) 00585 src_lons(4) = src_center_lon(n_add) 00586 00587 !*** 00588 !*** for consistency, we must make sure all lons are in 00589 !*** same 2pi interval 00590 !*** 00591 00592 vec1_lon = src_lons(1) - plon 00593 if (vec1_lon > pi) then 00594 src_lons(1) = src_lons(1) - pi2 00595 else if (vec1_lon < -pi) then 00596 src_lons(1) = src_lons(1) + pi2 00597 endif 00598 do n=2,4 00599 vec1_lon = src_lons(n) - src_lons(1) 00600 if (vec1_lon > pi) then 00601 src_lons(n) = src_lons(n) - pi2 00602 else if (vec1_lon < -pi) then 00603 src_lons(n) = src_lons(n) + pi2 00604 endif 00605 end do 00606 00607 corner_loop: do n=1,4 00608 next_n = MOD(n,4) + 1 00609 00610 !*** 00611 !*** here we take the cross product of the vector making 00612 !*** up each box side with the vector formed by the vertex 00613 !*** and search point. if all the cross products are 00614 !*** same sign, the point is contained in the box. 00615 !*** 00616 00617 vec1_lat = src_lats(next_n) - src_lats(n) 00618 vec1_lon = src_lons(next_n) - src_lons(n) 00619 vec2_lat = plat - src_lats(n) 00620 vec2_lon = plon - src_lons(n) 00621 00622 !*** 00623 !*** check for 0,2pi crossings 00624 !*** 00625 00626 if (vec1_lon > three*pih) then 00627 vec1_lon = vec1_lon - pi2 00628 else if (vec1_lon < -three*pih) then 00629 vec1_lon = vec1_lon + pi2 00630 endif 00631 if (vec2_lon > three*pih) then 00632 vec2_lon = vec2_lon - pi2 00633 else if (vec2_lon < -three*pih) then 00634 vec2_lon = vec2_lon + pi2 00635 endif 00636 00637 cross_product = vec1_lon*vec2_lat - vec2_lon*vec1_lat 00638 00639 !*** 00640 !*** if cross product is less than zero, this cell 00641 !*** doesn't work 00642 !*** 00643 00644 if (n == 1) cross_product_last = cross_product 00645 if (cross_product*cross_product_last < zero) then 00646 exit corner_loop 00647 else 00648 cross_product_last = cross_product 00649 endif 00650 00651 end do corner_loop 00652 00653 !*** 00654 !*** if cross products all positive, we found the location 00655 !*** 00656 00657 if (n > 4) then 00658 src_add(1) = srch_add 00659 src_add(2) = e_add 00660 src_add(3) = ne_add 00661 src_add(4) = n_add 00662 00663 ! Check if one point is masked; IF so, nearest-neighbour 00664 ! interpolation will be used 00665 00666 ntotmask = 0 00667 do ni=1,4 00668 if (.not. grid1_mask(src_add(ni)).and. 00669 & .not. lextrapdone) 00670 & ntotmask = ntotmask + 1 00671 end DO 00672 IF (ntotmask .gt. 0) src_add(1) = -src_add(1) 00673 00674 return 00675 endif 00676 00677 !*** 00678 !*** otherwise move on to next cell 00679 !*** 00680 00681 endif !bounding box check 00682 end do srch_loop 00683 00684 !*** 00685 !*** if no cell found, point is likely either in a box that straddles 00686 !*** either pole or is outside the grid. Put src_add = -1 so that 00687 !*** distance-weighted average of the 4 non-masked closest points 00688 !*** is done in calling routine. 00689 00690 src_add = -1 00691 00692 !----------------------------------------------------------------------- 00693 00694 end subroutine grid_search_bicub 00695 00696 !*********************************************************************** 00697 00698 subroutine store_link_bicub(dst_add, src_add, weights, nmap) 00699 00700 !----------------------------------------------------------------------- 00701 ! 00702 ! this routine stores the address and weight for four links 00703 ! associated with one destination point in the appropriate address 00704 ! and weight arrays and resizes those arrays if necessary. 00705 ! 00706 !----------------------------------------------------------------------- 00707 00708 !----------------------------------------------------------------------- 00709 ! 00710 ! input variables 00711 ! 00712 !----------------------------------------------------------------------- 00713 00714 integer (kind=int_kind), intent(in) :: 00715 & dst_add, ! address on destination grid 00716 & nmap ! identifies which direction for mapping 00717 00718 integer (kind=int_kind), dimension(4), intent(in) :: 00719 & src_add ! addresses on source grid 00720 00721 real (kind=dbl_kind), dimension(4,4), intent(in) :: 00722 & weights ! array of remapping weights for these links 00723 00724 !----------------------------------------------------------------------- 00725 ! 00726 ! local variables 00727 ! 00728 !----------------------------------------------------------------------- 00729 00730 integer (kind=int_kind) :: n, ! dummy index 00731 & num_links_old ! placeholder for old link number 00732 00733 !----------------------------------------------------------------------- 00734 ! 00735 ! increment number of links and check to see if remap arrays need 00736 ! to be increased to accomodate the new link. then store the 00737 ! link. 00738 ! 00739 !----------------------------------------------------------------------- 00740 00741 select case (nmap) 00742 case(1) 00743 00744 num_links_old = num_links_map1 00745 num_links_map1 = num_links_old + 4 00746 00747 if (num_links_map1 > max_links_map1) 00748 & call resize_remap_vars(1,resize_increment) 00749 00750 do n=1,4 00751 grid1_add_map1(num_links_old+n) = src_add(n) 00752 grid2_add_map1(num_links_old+n) = dst_add 00753 wts_map1 (:,num_links_old+n) = weights(:,n) 00754 end do 00755 00756 case(2) 00757 00758 num_links_old = num_links_map2 00759 num_links_map2 = num_links_old + 4 00760 00761 if (num_links_map2 > max_links_map2) 00762 & call resize_remap_vars(2,resize_increment) 00763 00764 do n=1,4 00765 grid1_add_map2(num_links_old+n) = dst_add 00766 grid2_add_map2(num_links_old+n) = src_add(n) 00767 wts_map2 (:,num_links_old+n) = weights(:,n) 00768 end do 00769 00770 end select 00771 00772 !----------------------------------------------------------------------- 00773 00774 end subroutine store_link_bicub 00775 00776 !*********************************************************************** 00777 00778 end module remap_bicubic 00779 00780 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1.7.4 
