|
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_bilinear module to calculate 00012 C bilinear remapping. 00013 C 00014 C** Interface: 00015 C --------- 00016 C *CALL* *remap_bilin** 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 ! bilinear interpolation. 00035 ! 00036 !----------------------------------------------------------------------- 00037 ! 00038 ! CVS:$Id: remap_bilinear.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_bilinear 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 implicit NONE 00074 00075 !----------------------------------------------------------------------- 00076 00077 integer (kind=int_kind), parameter :: 00078 & max_iter = 100 ! max iteration count for i,j iteration 00079 00080 real (kind=dbl_kind), parameter :: 00081 & converge = 1.e-10_dbl_kind ! convergence criterion 00082 00083 !*********************************************************************** 00084 00085 contains 00086 00087 !*********************************************************************** 00088 00089 subroutine remap_bilin (lextrapdone) 00090 00091 !----------------------------------------------------------------------- 00092 ! 00093 ! this routine computes the weights for a bilinear interpolation. 00094 ! 00095 !----------------------------------------------------------------------- 00096 00097 LOGICAL :: 00098 & 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 bilinear corners 00116 & src_lons, ! longitudes of four bilinear corners 00117 & wgts ! bilinear weights for four corners 00118 00119 real (kind=dbl_kind) :: 00120 & plat, plon, ! lat/lon coords of destination point 00121 & iguess, jguess, ! current guess for bilinear coordinate 00122 & deli, delj, ! corrections to i,j 00123 & dth1, dth2, dth3, ! some latitude differences 00124 & dph1, dph2, dph3, ! some longitude differences 00125 & dthp, dphp, ! difference between point and sw corner 00126 & mat1, mat2, mat3, mat4, ! matrix elements 00127 & determinant, ! matrix determinant 00128 & sum_wgts ! sum of weights for normalization 00129 00130 real (kind=dbl_kind) :: 00131 & coslat_dst, sinlat_dst, coslon_dst, sinlon_dst, 00132 & dist_min, distance, ! for computing dist-weighted avg 00133 & src_latsnn, arg 00134 00135 integer (kind=int_kind) :: min_add, max_add, srch_add, src_addnn 00136 00137 !----------------------------------------------------------------------- 00138 ! 00139 ! compute mappings from grid1 to grid2 00140 ! 00141 !----------------------------------------------------------------------- 00142 ! 00143 IF (nlogprt .GE. 2) THEN 00144 WRITE (UNIT = nulou,FMT = *)' ' 00145 WRITE (UNIT = nulou,FMT = *)'Entering routine remap_bilin' 00146 CALL FLUSH(nulou) 00147 ENDIF 00148 ! 00149 ll_nnei = .true. 00150 nmap = 1 00151 wgts=0 00152 src_add=0 00153 if (grid1_rank /= 2) then 00154 stop 'Can not do bilinear 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)) THEN 00164 cycle grid_loop1 00165 endif 00166 00167 plat = grid2_center_lat(dst_add) 00168 plon = grid2_center_lon(dst_add) 00169 00170 !*** 00171 !*** find nearest square of grid points on source grid 00172 !*** 00173 00174 call grid_search_bilin(src_add, src_lats, src_lons, 00175 & min_add, max_add, 00176 & plat, plon, grid1_dims, 00177 & grid1_center_lat, grid1_center_lon, 00178 & grid1_bound_box, bin_addr1, lextrapdone) 00179 00180 if (src_add(1) > 0) THEN 00181 00182 !*** 00183 !*** if the 4 surrounding points have been found and are 00184 !*** non-masked, do bilinear interpolation 00185 !*** 00186 00187 grid2_frac(dst_add) = one 00188 00189 00190 dth1 = src_lats(2) - src_lats(1) 00191 dth2 = src_lats(4) - src_lats(1) 00192 dth3 = src_lats(3) - src_lats(2) - dth2 00193 00194 dph1 = src_lons(2) - src_lons(1) 00195 dph2 = src_lons(4) - src_lons(1) 00196 dph3 = src_lons(3) - src_lons(2) 00197 00198 if (dph1 > three*pih) dph1 = dph1 - pi2 00199 if (dph2 > three*pih) dph2 = dph2 - pi2 00200 if (dph3 > three*pih) dph3 = dph3 - pi2 00201 if (dph1 < -three*pih) dph1 = dph1 + pi2 00202 if (dph2 < -three*pih) dph2 = dph2 + pi2 00203 if (dph3 < -three*pih) dph3 = dph3 + pi2 00204 00205 dph3 = dph3 - dph2 00206 00207 iguess = half 00208 jguess = half 00209 00210 iter_loop1: do iter=1,max_iter 00211 00212 dthp = plat - src_lats(1) - dth1*iguess - 00213 & dth2*jguess - dth3*iguess*jguess 00214 dphp = plon - src_lons(1) 00215 00216 if (dphp > three*pih) dphp = dphp - pi2 00217 if (dphp < -three*pih) dphp = dphp + pi2 00218 00219 dphp = dphp - dph1*iguess - dph2*jguess - 00220 & dph3*iguess*jguess 00221 00222 mat1 = dth1 + dth3*jguess 00223 mat2 = dth2 + dth3*iguess 00224 mat3 = dph1 + dph3*jguess 00225 mat4 = dph2 + dph3*iguess 00226 00227 determinant = mat1*mat4 - mat2*mat3 00228 00229 deli = (dthp*mat4 - mat2*dphp)/determinant 00230 delj = (mat1*dphp - dthp*mat3)/determinant 00231 00232 if (abs(deli) < converge .and. 00233 & abs(delj) < converge) exit iter_loop1 00234 00235 iguess = iguess + deli 00236 jguess = jguess + delj 00237 00238 end do iter_loop1 00239 00240 if (iter <= max_iter) then 00241 00242 !*** 00243 !*** successfully found i,j - compute weights 00244 !*** 00245 00246 wgts(1) = (one-iguess)*(one-jguess) 00247 wgts(2) = iguess*(one-jguess) 00248 wgts(3) = iguess*jguess 00249 wgts(4) = (one-iguess)*jguess 00250 00251 call store_link_bilin(dst_add, src_add, wgts, nmap) 00252 00253 ELSE 00254 write(nulou,*)'Point coords: ',plat,plon 00255 write(nulou,*)'Dest grid lats: ',src_lats 00256 write(nulou,*)'Dest grid lons: ',src_lons 00257 write(nulou,*)'Dest grid addresses: ',src_add 00258 write(nulou,*)'Current i,j : ',iguess, jguess 00259 write(nulou,*)'Iteration for i,j exceed max iteration count' 00260 stop 00261 endif 00262 00263 else if (src_add(1) < 0) THEN 00264 00265 !*** 00266 !*** Search for bilinear failed or at least one of the 4 00267 !*** neighbours was masked. Do distance-weighted average using 00268 !*** the non-masked points among the 4 closest ones. 00269 !*** 00270 00271 IF (nlogprt .eq. 2) then 00272 WRITE(nulou,*) ' ' 00273 WRITE(nulou,*) 00274 & 'WARNING: Cannot make bilinear interpolation for target point' 00275 & ,dst_add 00276 WRITE(nulou,*) 00277 & 'Using non-masked points among 4 nearest neighbors.' 00278 WRITE(nulou,*) ' ' 00279 ENDIF 00280 00281 !*** 00282 !*** Find the 4 closest points 00283 !*** 00284 00285 coslat_dst = cos(plat) 00286 sinlat_dst = sin(plat) 00287 coslon_dst = cos(plon) 00288 sinlon_dst = sin(plon) 00289 src_add = 0 00290 dist_min = bignum 00291 src_lats = bignum 00292 do srch_add = min_add,max_add 00293 arg = coslat_dst*cos(grid1_center_lat(srch_add))* 00294 & (coslon_dst*cos(grid1_center_lon(srch_add)) + 00295 & sinlon_dst*sin(grid1_center_lon(srch_add)))+ 00296 & sinlat_dst*sin(grid1_center_lat(srch_add)) 00297 IF (arg < -1.0d0) THEN 00298 arg = -1.0d0 00299 ELSE IF (arg > 1.0d0) THEN 00300 arg = 1.0d0 00301 END IF 00302 distance=acos(arg) 00303 00304 if (distance < dist_min) then 00305 sort_loop: do n=1,4 00306 if (distance < src_lats(n)) then 00307 do i=4,n+1,-1 00308 src_add (i) = src_add (i-1) 00309 src_lats(i) = src_lats(i-1) 00310 end do 00311 src_add (n) = srch_add 00312 src_lats(n) = distance 00313 dist_min = src_lats(4) 00314 exit sort_loop 00315 endif 00316 end do sort_loop 00317 endif 00318 end do 00319 00320 src_lons = one/(src_lats + tiny) 00321 distance = sum(src_lons) 00322 src_lats = src_lons/distance 00323 00324 !*** 00325 !*** Among 4 closest points, keep only the non-masked ones 00326 !*** 00327 00328 icount = 0 00329 do n=1,4 00330 if (grid1_mask(src_add(n)) .or. 00331 & (.not. grid1_mask(src_add(n)) .and. lextrapdone)) then 00332 icount = icount + 1 00333 else 00334 src_lats(n) = zero 00335 endif 00336 end do 00337 if (icount > 0) then 00338 !*** renormalize weights 00339 sum_wgts = sum(src_lats) 00340 wgts(1) = src_lats(1)/sum_wgts 00341 wgts(2) = src_lats(2)/sum_wgts 00342 wgts(3) = src_lats(3)/sum_wgts 00343 wgts(4) = src_lats(4)/sum_wgts 00344 00345 grid2_frac(dst_add) = one 00346 call store_link_bilin(dst_add, src_add, wgts, nmap) 00347 ELSE 00348 IF (ll_nnei .eqv. .true. ) then 00349 IF (nlogprt .ge. 2) THEN 00350 WRITE(nulou,*) ' ' 00351 WRITE(nulou,*) 00352 & 'All 4 surrounding source grid points are masked' 00353 WRITE(nulou,*) 'for target point ',dst_add 00354 WRITE(nulou,*) 'with longitude and latitude', 00355 & plon, plat 00356 WRITE(nulou,*) 00357 & 'Using the nearest non-masked neighbour.' 00358 CALL FLUSH(nulou) 00359 ENDIF 00360 C 00361 src_latsnn = bignum 00362 !cdir novector 00363 do srch_add = min_add,max_add 00364 if (grid1_mask(srch_add) .or. 00365 & (.not. grid1_mask(srch_add) .and. lextrapdone)) THEN 00366 arg = coslat_dst*cos(grid1_center_lat(srch_add))* 00367 & (coslon_dst*cos(grid1_center_lon(srch_add)) + 00368 & sinlon_dst*sin(grid1_center_lon(srch_add)))+ 00369 & sinlat_dst*sin(grid1_center_lat(srch_add)) 00370 IF (arg < -1.0d0) THEN 00371 arg = -1.0d0 00372 ELSE IF (arg > 1.0d0) THEN 00373 arg = 1.0d0 00374 END IF 00375 distance=acos(arg) 00376 00377 if (distance < src_latsnn) then 00378 src_addnn = srch_add 00379 src_latsnn = distance 00380 endif 00381 endif 00382 end DO 00383 IF (nlogprt .ge. 2) THEN 00384 WRITE(nulou,*) ' ' 00385 WRITE(nulou,*) 00386 & 'Nearest non masked neighbour is source point ' 00387 & ,src_addnn 00388 WRITE(nulou,*) 'with longitude and latitude', 00389 & grid1_center_lon(src_addnn), 00390 & grid1_center_lat(src_addnn) 00391 ENDIF 00392 wgts(1) = 1. 00393 wgts(2) = 0. 00394 wgts(3) = 0. 00395 wgts(4) = 0. 00396 src_add(1) = src_addnn 00397 src_add(2) = 0 00398 src_add(3) = 0 00399 src_add(4) = 0 00400 00401 grid2_frac(dst_add) = one 00402 call store_link_bilin(dst_add, src_add, wgts, nmap) 00403 00404 ENDIF 00405 ENDIF 00406 ENDIF 00407 end do grid_loop1 00408 ! 00409 IF (nlogprt .GE. 2) THEN 00410 WRITE (UNIT = nulou,FMT = *)' ' 00411 WRITE (UNIT = nulou,FMT = *)'Leaving routine remap_bilin' 00412 CALL FLUSH(nulou) 00413 ENDIF 00414 ! 00415 end subroutine remap_bilin 00416 00417 !*********************************************************************** 00418 00419 subroutine grid_search_bilin(src_add, src_lats, src_lons, 00420 & min_add, max_add, 00421 & plat, plon, src_grid_dims, 00422 & src_center_lat, src_center_lon, 00423 & src_grid_bound_box, 00424 & src_bin_add, lextrapdone) 00425 00426 !----------------------------------------------------------------------- 00427 ! 00428 ! this routine finds the location of the search point plat, plon 00429 ! in the source grid and returns the corners needed for a bilinear 00430 ! interpolation. 00431 ! 00432 !----------------------------------------------------------------------- 00433 00434 !----------------------------------------------------------------------- 00435 ! 00436 ! output variables 00437 ! 00438 !----------------------------------------------------------------------- 00439 00440 integer (kind=int_kind), dimension(4), intent(out) :: 00441 & src_add ! address of each corner point enclosing P 00442 00443 real (kind=dbl_kind), dimension(4), intent(out) :: 00444 & src_lats, ! latitudes of the four corner points 00445 & src_lons ! longitudes of the four corner points 00446 00447 integer (kind=int_kind) :: min_add, max_add 00448 00449 !----------------------------------------------------------------------- 00450 ! 00451 ! input variables 00452 ! 00453 !----------------------------------------------------------------------- 00454 00455 real (kind=dbl_kind), intent(in) :: 00456 & plat, ! latitude of the search point 00457 & plon ! longitude of the search point 00458 00459 integer (kind=int_kind), dimension(2), intent(in) :: 00460 & src_grid_dims ! size of each src grid dimension 00461 00462 real (kind=dbl_kind), dimension(:), intent(in) :: 00463 & src_center_lat, ! latitude of each src grid center 00464 & src_center_lon ! longitude of each src grid center 00465 00466 real (kind=dbl_kind), dimension(:,:), intent(in) :: 00467 & src_grid_bound_box ! bound box for source grid 00468 00469 integer (kind=int_kind), dimension(:,:), intent(in) :: 00470 & src_bin_add ! latitude bins for restricting 00471 00472 LOGICAL :: 00473 & lextrapdone ! logical, true if EXTRAP done on field 00474 !----------------------------------------------------------------------- 00475 ! 00476 ! local variables 00477 ! 00478 !----------------------------------------------------------------------- 00479 00480 integer (kind=int_kind) :: n, next_n, srch_add, ni, ! dummy indices 00481 & nx, ny, ntotmask, ! dimensions of src grid 00482 & i, j, jp1, ip1, n_add, e_add, ne_add ! addresses 00483 00484 real (kind=dbl_kind) :: ! vectors for cross-product check 00485 & vec1_lat, vec1_lon, 00486 & vec2_lat, vec2_lon, cross_product, cross_product_last 00487 00488 !----------------------------------------------------------------------- 00489 ! 00490 ! restrict search first using bins 00491 ! 00492 !----------------------------------------------------------------------- 00493 00494 src_add = 0 00495 00496 min_add = size(src_center_lat) 00497 max_add = 1 00498 do n=1,num_srch_bins 00499 if (plat >= bin_lats(1,n) .and. plat <= bin_lats(2,n) .and. 00500 & plon >= bin_lons(1,n) .and. plon <= bin_lons(2,n)) then 00501 min_add = min(min_add, src_bin_add(1,n)) 00502 max_add = max(max_add, src_bin_add(2,n)) 00503 endif 00504 end do 00505 00506 !----------------------------------------------------------------------- 00507 ! 00508 ! now perform a more detailed search 00509 ! 00510 !----------------------------------------------------------------------- 00511 00512 nx = src_grid_dims(1) 00513 ny = src_grid_dims(2) 00514 00515 srch_loop: do srch_add = min_add,max_add 00516 00517 !*** first check bounding box 00518 00519 if (plat <= src_grid_bound_box(2,srch_add) .and. 00520 & plat >= src_grid_bound_box(1,srch_add) .and. 00521 & plon <= src_grid_bound_box(4,srch_add) .and. 00522 & plon >= src_grid_bound_box(3,srch_add)) then 00523 00524 !*** 00525 !*** we are within bounding box so get really serious 00526 !*** 00527 00528 !*** determine neighbor addresses 00529 00530 j = (srch_add - 1)/nx +1 00531 i = srch_add - (j-1)*nx 00532 00533 !*** find ip1 00534 !*** Note: I do not want to take into account the number 00535 !*** of overlapping grid points, as the underlying cell 00536 !*** will be found in all cases if the grid overlaps. 00537 00538 if (i < nx) then 00539 ip1 = i + 1 00540 else 00541 ip1 = 1 00542 endif 00543 00544 if (j < ny) then 00545 jp1 = j+1 00546 else 00547 jp1 = 1 00548 endif 00549 00550 n_add = (jp1 - 1)*nx + i 00551 e_add = (j - 1)*nx + ip1 00552 ne_add = (jp1 - 1)*nx + ip1 00553 00554 src_lats(1) = src_center_lat(srch_add) 00555 src_lats(2) = src_center_lat(e_add) 00556 src_lats(3) = src_center_lat(ne_add) 00557 src_lats(4) = src_center_lat(n_add) 00558 00559 src_lons(1) = src_center_lon(srch_add) 00560 src_lons(2) = src_center_lon(e_add) 00561 src_lons(3) = src_center_lon(ne_add) 00562 src_lons(4) = src_center_lon(n_add) 00563 00564 !*** 00565 !*** for consistency, we must make sure all lons are in 00566 !*** same 2pi interval 00567 !*** 00568 00569 vec1_lon = src_lons(1) - plon 00570 if (vec1_lon > pi) then 00571 src_lons(1) = src_lons(1) - pi2 00572 else if (vec1_lon < -pi) then 00573 src_lons(1) = src_lons(1) + pi2 00574 endif 00575 do n=2,4 00576 vec1_lon = src_lons(n) - src_lons(1) 00577 if (vec1_lon > pi) then 00578 src_lons(n) = src_lons(n) - pi2 00579 else if (vec1_lon < -pi) then 00580 src_lons(n) = src_lons(n) + pi2 00581 endif 00582 end do 00583 00584 corner_loop: do n=1,4 00585 next_n = MOD(n,4) + 1 00586 00587 !*** 00588 !*** here we take the cross product of the vector making 00589 !*** up each box side with the vector formed by the vertex 00590 !*** and search point. if all the cross products are 00591 !*** positive, the point is contained in the box. 00592 !*** 00593 00594 vec1_lat = src_lats(next_n) - src_lats(n) 00595 vec1_lon = src_lons(next_n) - src_lons(n) 00596 vec2_lat = plat - src_lats(n) 00597 vec2_lon = plon - src_lons(n) 00598 00599 !*** 00600 !*** check for 0,2pi crossings 00601 !*** 00602 00603 if (vec1_lon > three*pih) then 00604 vec1_lon = vec1_lon - pi2 00605 else if (vec1_lon < -three*pih) then 00606 vec1_lon = vec1_lon + pi2 00607 endif 00608 if (vec2_lon > three*pih) then 00609 vec2_lon = vec2_lon - pi2 00610 else if (vec2_lon < -three*pih) then 00611 vec2_lon = vec2_lon + pi2 00612 endif 00613 00614 cross_product = vec1_lon*vec2_lat - vec2_lon*vec1_lat 00615 00616 !*** 00617 !*** if cross product is less than zero, this cell 00618 !*** doesn't work 00619 !*** 00620 00621 if (n == 1) cross_product_last = cross_product 00622 if (cross_product*cross_product_last < zero) 00623 & exit corner_loop 00624 cross_product_last = cross_product 00625 00626 end do corner_loop 00627 00628 !*** 00629 !*** if cross products all same sign, we found the location 00630 !*** 00631 00632 if (n > 4) then 00633 src_add(1) = srch_add 00634 src_add(2) = e_add 00635 src_add(3) = ne_add 00636 src_add(4) = n_add 00637 00638 ! Check if one point is masked; IF so, nearest-neighbour 00639 ! interpolation will be used 00640 00641 ntotmask = 0 00642 do ni=1,4 00643 if (.not. grid1_mask(src_add(ni)).and. 00644 & .not. lextrapdone) 00645 & ntotmask = ntotmask + 1 00646 end DO 00647 IF (ntotmask .gt. 0) src_add(1) = -src_add(1) 00648 00649 return 00650 endif 00651 00652 !*** 00653 !*** otherwise move on to next cell 00654 !*** 00655 00656 endif !bounding box check 00657 end do srch_loop 00658 00659 !*** 00660 !*** if no cell found, point is likely either in a box that straddles 00661 !*** either pole or is outside the grid. Put src_add = -1 so that 00662 !*** distance-weighted average of the 4 non-masked closest points 00663 !*** is done in calling routine. 00664 00665 src_add = -1 00666 00667 !----------------------------------------------------------------------- 00668 00669 end subroutine grid_search_bilin 00670 00671 !*********************************************************************** 00672 00673 subroutine store_link_bilin(dst_add, src_add, weights, nmap) 00674 00675 !----------------------------------------------------------------------- 00676 ! 00677 ! this routine stores the address and weight for four links 00678 ! associated with one destination point in the appropriate address 00679 ! and weight arrays and resizes those arrays if necessary. 00680 ! 00681 !----------------------------------------------------------------------- 00682 00683 !----------------------------------------------------------------------- 00684 ! 00685 ! input variables 00686 ! 00687 !----------------------------------------------------------------------- 00688 00689 integer (kind=int_kind), intent(in) :: 00690 & dst_add, ! address on destination grid 00691 & nmap ! identifies which direction for mapping 00692 00693 integer (kind=int_kind), dimension(4), intent(in) :: 00694 & src_add ! addresses on source grid 00695 00696 real (kind=dbl_kind), dimension(4), intent(in) :: 00697 & weights ! array of remapping weights for these links 00698 00699 !----------------------------------------------------------------------- 00700 ! 00701 ! local variables 00702 ! 00703 !----------------------------------------------------------------------- 00704 00705 integer (kind=int_kind) :: n, ! dummy index 00706 & num_links_old ! placeholder for old link number 00707 00708 !----------------------------------------------------------------------- 00709 ! 00710 ! increment number of links and check to see if remap arrays need 00711 ! to be increased to accomodate the new link. then store the 00712 ! link. 00713 ! 00714 !----------------------------------------------------------------------- 00715 00716 select case (nmap) 00717 case(1) 00718 00719 num_links_old = num_links_map1 00720 num_links_map1 = num_links_old + 4 00721 00722 if (num_links_map1 > max_links_map1) 00723 & call resize_remap_vars(1,resize_increment) 00724 00725 do n=1,4 00726 grid1_add_map1(num_links_old+n) = src_add(n) 00727 grid2_add_map1(num_links_old+n) = dst_add 00728 wts_map1 (1,num_links_old+n) = weights(n) 00729 end do 00730 00731 case(2) 00732 00733 num_links_old = num_links_map2 00734 num_links_map2 = num_links_old + 4 00735 00736 if (num_links_map2 > max_links_map2) 00737 & call resize_remap_vars(2,resize_increment) 00738 00739 do n=1,4 00740 grid1_add_map2(num_links_old+n) = dst_add 00741 grid2_add_map2(num_links_old+n) = src_add(n) 00742 wts_map2 (1,num_links_old+n) = weights(n) 00743 end do 00744 00745 end select 00746 00747 !----------------------------------------------------------------------- 00748 00749 end subroutine store_link_bilin 00750 00751 !*********************************************************************** 00752 00753 end module remap_bilinear 00754 00755 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1.7.4 
