- Timestamp:
- 2020-06-03T16:26:23+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
utils/tools_dev_r12970_AGRIF_CMEMS/DOMAINcfg/src/bilinear_interp.f90
r12414 r13024 1 ! 1 2 MODULE bilinear_interp 2 3 END MODULE 3 ! 4 USE agrif_modutil 5 ! 6 !************************************************************************ 7 ! * 8 ! MODULE BILINEAR INTERP * 9 ! * 10 ! bilinear interpolation routines from SCRIP package * 11 ! * 12 !http://climate.lanl.gov/Software/SCRIP/ * 13 ! * 14 !Bilinear remapping * 15 ! * 16 !************************************************************************ 17 ! 18 !----------------------------------------------------------------------- 19 IMPLICIT NONE 20 21 !----------------------------------------------------------------------- 22 ! variables that describe each grid 23 !----------------------------------------------------------------------- 24 ! 25 INTEGER :: grid1_size,grid2_size,grid1_rank, grid2_rank 26 ! 27 INTEGER, DIMENSION(:), ALLOCATABLE :: grid1_dims, grid2_dims 28 ! 29 30 !----------------------------------------------------------------------- 31 ! grid coordinates and masks 32 !----------------------------------------------------------------------- 33 ! 34 LOGICAL, DIMENSION(:), ALLOCATABLE :: grid1_mask,grid2_mask 35 ! each grid center in radians 36 REAL*8,DIMENSION(:), ALLOCATABLE :: & 37 grid1_center_lat, & 38 grid1_center_lon, & 39 grid2_center_lat, & 40 grid2_center_lon, & 41 grid1_frac, & ! fractional area of grid cells 42 grid2_frac ! participating in remapping 43 ! 44 ! lat/lon bounding box for use in restricting grid searches 45 ! 46 REAL*8,DIMENSION(:,:), ALLOCATABLE :: grid1_bound_box,grid2_bound_box 47 ! 48 !----------------------------------------------------------------------- 49 ! bins for restricting searches 50 !----------------------------------------------------------------------- 51 ! 52 ! num of bins for restricted srch 53 INTEGER, PARAMETER :: num_srch_bins = 90 54 ! 55 ! min,max adds for grid cells in this lat bin 56 ! 57 INTEGER,DIMENSION(:,:), ALLOCATABLE :: bin_addr1,bin_addr2 58 ! 59 ! min,max longitude for each search bin 60 ! 61 REAL*8, DIMENSION(:,:), ALLOCATABLE :: bin_lats,bin_lons 62 63 REAL*8, PARAMETER :: zero = 0.0, & 64 one = 1.0, & 65 two = 2.0, & 66 three = 3.0, & 67 four = 4.0, & 68 five = 5.0, & 69 half = 0.5, & 70 quart = 0.25, & 71 bignum = 1.e+20, & 72 tiny = 1.e-14, & 73 pi = 3.14159265359, & 74 pi2 = two*pi, & 75 pih = half*pi 76 77 REAL*8, PARAMETER :: deg2rad = pi/180. 78 ! 79 ! max iteration count for i,j iteration 80 ! 81 INTEGER , PARAMETER :: max_iter = 100 82 ! 83 ! convergence criterion 84 ! 85 REAL*8, PARAMETER :: converge = 1.e-10 86 87 88 ! 89 INTEGER, PARAMETER :: norm_opt_none = 1 & 90 ,norm_opt_dstarea = 2 & 91 ,norm_opt_frcarea = 3 92 ! 93 INTEGER, PARAMETER :: map_type_conserv = 1 & 94 ,map_type_bilinear = 2 & 95 ,map_type_bicubic = 3 & 96 ,map_type_distwgt = 4 97 ! 98 INTEGER :: max_links_map1 & ! current size of link arrays 99 ,num_links_map1 & ! actual number of links for remapping 100 ,max_links_map2 & ! current size of link arrays 101 ,num_links_map2 & ! actual number of links for remapping 102 ,num_maps & ! num of remappings for this grid pair 103 ,num_wts & ! num of weights used in remapping 104 ,map_type & ! identifier for remapping method 105 ,norm_opt & ! option for normalization (conserv only) 106 ,resize_increment ! default amount to increase array size 107 108 INTEGER , DIMENSION(:), ALLOCATABLE :: & 109 grid1_add_map1, & ! grid1 address for each link in mapping 1 110 grid2_add_map1, & ! grid2 address for each link in mapping 1 111 grid1_add_map2, & ! grid1 address for each link in mapping 2 112 grid2_add_map2 ! grid2 address for each link in mapping 2 113 114 REAL*8, DIMENSION(:,:), ALLOCATABLE :: & 115 wts_map1, & ! map weights for each link (num_wts,max_links) 116 wts_map2 ! map weights for each link (num_wts,max_links) 117 ! 118 CONTAINS 119 ! 120 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 121 !************************************************************************ 122 ! SUBROUTINE GRID_INIT 123 !************************************************************************ 124 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 125 ! 126 SUBROUTINE get_remap_matrix(grid1_lat,grid2_lat,grid1_lon,grid2_lon,mask, & 127 remap_matrix,source_add,destination_add) 128 ! 129 !----------------------------------------------------------------------- 130 !this routine makes any necessary changes (e.g. for 0,2pi longitude range) 131 !----------------------------------------------------------------------- 132 ! 133 REAL*8,DIMENSION(:,:) :: grid1_lat,grid2_lat,grid1_lon,grid2_lon 134 LOGICAL,DIMENSION(:,:) :: mask 135 ! 136 INTEGER,DIMENSION(:),ALLOCATABLE :: source_add,destination_add 137 REAL*8,DIMENSION(:,:), ALLOCATABLE :: remap_matrix 138 ! 139 !----------------------------------------------------------------------- 140 ! local variables 141 !----------------------------------------------------------------------- 142 ! 143 INTEGER :: n,nele,i,j,ip1,jp1,n_add,e_add,ne_add,nx,ny 144 INTEGER :: xpos,ypos 145 ! 146 ! integer mask 147 ! 148 INTEGER, DIMENSION(:), ALLOCATABLE :: imask 149 ! 150 ! lat/lon intervals for search bins 151 ! 152 REAL*8 :: dlat,dlon 153 ! 154 ! temps for computing bounding boxes 155 ! 156 REAL*8, DIMENSION(4) :: tmp_lats, tmp_lons 157 ! 158 ! write(*,*)'proceed to Bilinear interpolation ...' 159 ! 160 ! IF(ALLOCATED(wts_map1)) DEALLOCATE(wts_map1) 161 ! IF(ALLOCATED(grid1_add_map1)) DEALLOCATE(grid1_add_map1) 162 ! IF(ALLOCATED(grid2_add_map1)) DEALLOCATE(grid2_add_map1) 163 164 165 ! 166 ALLOCATE(grid1_dims(2),grid2_dims(2)) 167 ! 168 grid1_dims(1) = SIZE(grid1_lat,2) 169 grid1_dims(2) = SIZE(grid1_lat,1) 170 grid2_dims(1) = SIZE(grid2_lat,2) 171 grid2_dims(2) = SIZE(grid2_lat,1) 172 grid1_size = SIZE(grid1_lat,2) * SIZE(grid1_lat,1) 173 grid2_size = SIZE(grid2_lat,2) * SIZE(grid2_lat,1) 174 ! 175 !----------------------------------------------------------------------- 176 ! allocate grid coordinates/masks and read data 177 !----------------------------------------------------------------------- 178 ! 179 ALLOCATE( grid2_mask(grid2_size), & 180 grid1_bound_box (4,grid1_size), & 181 grid2_bound_box (4,grid2_size), & 182 grid1_frac (grid1_size), & 183 grid2_frac (grid2_size)) 184 ALLOCATE(imask(grid1_size)) 185 ! 186 ! 187 grid1_frac = zero 188 grid2_frac = zero 189 190 ! 191 ! 2D array -> 1D array 192 ! 193 ALLOCATE(grid1_center_lat(SIZE(grid1_lat,1)*SIZE(grid1_lat,2))) 194 CALL tab2Dto1D(grid1_lat,grid1_center_lat) 195 196 ALLOCATE(grid1_center_lon(SIZE(grid1_lon,1)*SIZE(grid1_lon,2))) 197 CALL tab2Dto1D(grid1_lon,grid1_center_lon) 198 199 ALLOCATE(grid2_center_lat(SIZE(grid2_lat,1)*SIZE(grid2_lat,2))) 200 CALL tab2Dto1D(grid2_lat,grid2_center_lat) 201 202 ALLOCATE(grid2_center_lon(SIZE(grid2_lon,1)*SIZE(grid2_lon,2))) 203 CALL tab2Dto1D(grid2_lon,grid2_center_lon) 204 205 ALLOCATE(grid1_mask(SIZE(grid1_lat,1)*SIZE(grid1_lat,2))) 206 CALL logtab2Dto1D(mask,grid1_mask) 207 ! 208 ! Write(*,*) ,'grid1_mask = ',grid1_mask 209 ! 210 ! degrees to radian 211 ! 212 grid1_center_lat = grid1_center_lat*deg2rad 213 grid1_center_lon = grid1_center_lon*deg2rad 214 grid2_center_lat = grid2_center_lat*deg2rad 215 grid2_center_lon = grid2_center_lon*deg2rad 216 217 !----------------------------------------------------------------------- 218 ! convert longitudes to 0,2pi interval 219 !----------------------------------------------------------------------- 220 221 WHERE (grid1_center_lon .GT. pi2) grid1_center_lon = & 222 grid1_center_lon - pi2 223 WHERE (grid1_center_lon .LT. zero) grid1_center_lon = & 224 grid1_center_lon + pi2 225 WHERE (grid2_center_lon .GT. pi2) grid2_center_lon = & 226 grid2_center_lon - pi2 227 WHERE (grid2_center_lon .LT. zero) grid2_center_lon = & 228 grid2_center_lon + pi2 229 230 !----------------------------------------------------------------------- 231 ! 232 ! make sure input latitude range is within the machine values 233 ! for +/- pi/2 234 ! 235 !----------------------------------------------------------------------- 236 237 WHERE (grid1_center_lat > pih) grid1_center_lat = pih 238 WHERE (grid1_center_lat < -pih) grid1_center_lat = -pih 239 WHERE (grid2_center_lat > pih) grid2_center_lat = pih 240 WHERE (grid2_center_lat < -pih) grid2_center_lat = -pih 241 242 !----------------------------------------------------------------------- 243 ! 244 ! compute bounding boxes for restricting future grid searches 245 ! 246 !----------------------------------------------------------------------- 247 ! 248 nx = grid1_dims(1) 249 ny = grid1_dims(2) 250 251 DO n=1,grid1_size 252 253 !*** find N,S and NE points to this grid point 254 255 j = (n - 1)/nx +1 256 i = n - (j-1)*nx 257 258 IF (i < nx) THEN 259 ip1 = i + 1 260 ELSE 261 !*** assume cyclic 262 ip1 = 1 263 !*** but if it is not, correct 264 e_add = (j - 1)*nx + ip1 265 IF (ABS(grid1_center_lat(e_add) - & 266 grid1_center_lat(n )) > pih) THEN 267 ip1 = i 268 ENDIF 269 ip1=nx 270 ENDIF 271 272 IF (j < ny) THEN 273 jp1 = j+1 274 ELSE 275 !*** assume cyclic 276 jp1 = 1 277 !*** but if it is not, correct 278 n_add = (jp1 - 1)*nx + i 279 IF (ABS(grid1_center_lat(n_add) - & 280 grid1_center_lat(n )) > pih) THEN 281 jp1 = j 282 ENDIF 283 jp1=ny 284 ENDIF 285 286 n_add = (jp1 - 1)*nx + i 287 e_add = (j - 1)*nx + ip1 288 ne_add = (jp1 - 1)*nx + ip1 289 290 !*** find N,S and NE lat/lon coords and check bounding box 291 292 tmp_lats(1) = grid1_center_lat(n) 293 tmp_lats(2) = grid1_center_lat(e_add) 294 tmp_lats(3) = grid1_center_lat(ne_add) 295 tmp_lats(4) = grid1_center_lat(n_add) 296 297 tmp_lons(1) = grid1_center_lon(n) 298 tmp_lons(2) = grid1_center_lon(e_add) 299 tmp_lons(3) = grid1_center_lon(ne_add) 300 tmp_lons(4) = grid1_center_lon(n_add) 301 302 grid1_bound_box(1,n) = MINVAL(tmp_lats) 303 grid1_bound_box(2,n) = MAXVAL(tmp_lats) 304 305 grid1_bound_box(3,n) = MINVAL(tmp_lons) 306 grid1_bound_box(4,n) = MAXVAL(tmp_lons) 307 END DO 308 309 nx = grid2_dims(1) 310 ny = grid2_dims(2) 311 312 DO n=1,grid2_size 313 314 !*** find N,S and NE points to this grid point 315 316 j = (n - 1)/nx +1 317 i = n - (j-1)*nx 318 319 IF (i < nx) THEN 320 ip1 = i + 1 321 ELSE 322 !*** assume cyclic 323 ip1 = 1 324 !*** but if it is not, correct 325 e_add = (j - 1)*nx + ip1 326 IF (ABS(grid2_center_lat(e_add) - & 327 grid2_center_lat(n )) > pih) THEN 328 ip1 = i 329 ENDIF 330 ENDIF 331 332 IF (j < ny) THEN 333 jp1 = j+1 334 ELSE 335 !*** assume cyclic 336 jp1 = 1 337 !*** but if it is not, correct 338 n_add = (jp1 - 1)*nx + i 339 IF (ABS(grid2_center_lat(n_add) - & 340 grid2_center_lat(n )) > pih) THEN 341 jp1 = j 342 ENDIF 343 ENDIF 344 345 n_add = (jp1 - 1)*nx + i 346 e_add = (j - 1)*nx + ip1 347 ne_add = (jp1 - 1)*nx + ip1 348 349 !*** find N,S and NE lat/lon coords and check bounding box 350 351 tmp_lats(1) = grid2_center_lat(n) 352 tmp_lats(2) = grid2_center_lat(e_add) 353 tmp_lats(3) = grid2_center_lat(ne_add) 354 tmp_lats(4) = grid2_center_lat(n_add) 355 356 tmp_lons(1) = grid2_center_lon(n) 357 tmp_lons(2) = grid2_center_lon(e_add) 358 tmp_lons(3) = grid2_center_lon(ne_add) 359 tmp_lons(4) = grid2_center_lon(n_add) 360 361 grid2_bound_box(1,n) = MINVAL(tmp_lats) 362 grid2_bound_box(2,n) = MAXVAL(tmp_lats) 363 grid2_bound_box(3,n) = MINVAL(tmp_lons) 364 grid2_bound_box(4,n) = MAXVAL(tmp_lons) 365 END DO 366 ! 367 ! 368 ! 369 WHERE (ABS(grid1_bound_box(4,:) - grid1_bound_box(3,:)) > pi) 370 grid1_bound_box(3,:) = zero 371 grid1_bound_box(4,:) = pi2 372 END WHERE 373 374 WHERE (ABS(grid2_bound_box(4,:) - grid2_bound_box(3,:)) > pi) 375 grid2_bound_box(3,:) = zero 376 grid2_bound_box(4,:) = pi2 377 END WHERE 378 379 !*** 380 !*** try to check for cells that overlap poles 381 !*** 382 383 WHERE (grid1_center_lat > grid1_bound_box(2,:)) & 384 grid1_bound_box(2,:) = pih 385 386 WHERE (grid1_center_lat < grid1_bound_box(1,:)) & 387 grid1_bound_box(1,:) = -pih 388 389 WHERE (grid2_center_lat > grid2_bound_box(2,:)) & 390 grid2_bound_box(2,:) = pih 391 392 WHERE (grid2_center_lat < grid2_bound_box(1,:)) & 393 grid2_bound_box(1,:) = -pih 394 395 !----------------------------------------------------------------------- 396 ! set up and assign address ranges to search bins in order to 397 ! further restrict later searches 398 !----------------------------------------------------------------------- 399 400 ALLOCATE(bin_addr1(2,num_srch_bins)) 401 ALLOCATE(bin_addr2(2,num_srch_bins)) 402 ALLOCATE(bin_lats (2,num_srch_bins)) 403 ALLOCATE(bin_lons (2,num_srch_bins)) 404 405 dlat = pi/num_srch_bins 406 407 DO n=1,num_srch_bins 408 bin_lats(1,n) = (n-1)*dlat - pih 409 bin_lats(2,n) = n*dlat - pih 410 bin_lons(1,n) = zero 411 bin_lons(2,n) = pi2 412 bin_addr1(1,n) = grid1_size + 1 413 bin_addr1(2,n) = 0 414 bin_addr2(1,n) = grid2_size + 1 415 bin_addr2(2,n) = 0 416 END DO 417 418 DO nele=1,grid1_size 419 DO n=1,num_srch_bins 420 IF (grid1_bound_box(1,nele) <= bin_lats(2,n) .AND. & 421 grid1_bound_box(2,nele) >= bin_lats(1,n)) THEN 422 bin_addr1(1,n) = MIN(nele,bin_addr1(1,n)) 423 bin_addr1(2,n) = MAX(nele,bin_addr1(2,n)) 424 ENDIF 425 END DO 426 END DO 427 428 DO nele=1,grid2_size 429 DO n=1,num_srch_bins 430 IF (grid2_bound_box(1,nele) <= bin_lats(2,n) .AND. & 431 grid2_bound_box(2,nele) >= bin_lats(1,n)) THEN 432 bin_addr2(1,n) = MIN(nele,bin_addr2(1,n)) 433 bin_addr2(2,n) = MAX(nele,bin_addr2(2,n)) 434 ENDIF 435 END DO 436 END DO 437 ! 438 CALL init_remap_vars 439 CALL remap_bilin 440 441 ALLOCATE(remap_matrix(SIZE(wts_map1,1),SIZE(wts_map1,2)), & 442 source_add(SIZE(grid1_add_map1)), & 443 destination_add(SIZE(grid2_add_map1))) 444 445 DO j = 1,SIZE(wts_map1,2) 446 DO i = 1,SIZE(wts_map1,1) 447 448 remap_matrix(i,j) = wts_map1(i,j) 449 450 END DO 451 END DO 452 453 454 source_add(:) = grid1_add_map1(:) 455 destination_add(:) = grid2_add_map1(:) 456 ! 457 WHERE(destination_add == 0) 458 destination_add = 1 459 END WHERE 460 461 WHERE(source_add == 0) 462 source_add = 1 463 END WHERE 464 ! 465 !DEALLOCATE(grid1_bound_box,grid2_bound_box,grid1_center_lat,grid1_center_lon) 466 !DEALLOCATE(grid2_center_lat,grid2_center_lon,grid2_add_map1,grid1_add_map1,wts_map1) 467 !DEALLOCATE(grid1_frac,grid2_frac,grid1_dims,grid2_dims,grid2_mask,imask) 468 !DEALLOCATE(bin_addr1,bin_addr2,bin_lats,bin_lons) 469 !DEALLOCATE(grid1_mask) 470 ! 471 !----------------------------------------------------------------------- 472 473 END SUBROUTINE get_remap_matrix 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 !*********************************************************************** 499 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 500 !************************************************************************ 501 ! SUBROUTINE REMAP_BILINEAR 502 !************************************************************************ 503 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 504 505 SUBROUTINE remap_bilin 506 507 !----------------------------------------------------------------------- 508 ! this routine computes the weights for a bilinear interpolation. 509 !----------------------------------------------------------------------- 510 511 !----------------------------------------------------------------------- 512 ! local variables 513 !----------------------------------------------------------------------- 514 515 INTEGER :: n,icount,dst_add,iter,nmap,nbmasked 516 ! 517 ! address for the four source points 518 ! 519 INTEGER, DIMENSION(4) :: src_add,involved_pts 520 INTEGER, DIMENSION(1) :: minlon 521 INTEGER, DIMENSION(1) :: minlat 522 REAL*8, DIMENSION(4) :: distx,disty 523 REAL*8 :: normalize 524 ! 525 ! latitudes longitudes of four bilinear corners 526 ! 527 REAL*8, DIMENSION(4) :: src_lats,src_lons 528 ! 529 ! bilinear weights for four corners 530 ! 531 REAL*8, DIMENSION(4) :: wgts 532 ! 533 REAL*8 :: & 534 plat, plon, & ! lat/lon coords of destination point 535 iguess, jguess, & ! current guess for bilinear coordinate 536 thguess, phguess, & ! current guess for lat/lon coordinate 537 deli, delj, & ! corrections to i,j 538 dth1, dth2, dth3, & ! some latitude differences 539 dph1, dph2, dph3, & ! some longitude differences 540 dthp, dphp, & ! difference between point and sw corner 541 mat1, mat2, mat3, mat4, & ! matrix elements 542 determinant, & ! matrix determinant 543 sum_wgts ! sum of weights for normalization 544 545 INTEGER lastsrc_add 546 REAL*8,DIMENSION(:),ALLOCATABLE :: cos_grid1_center_lat 547 REAL*8,DIMENSION(:),ALLOCATABLE :: cos_grid1_center_lon 548 REAL*8,DIMENSION(:),ALLOCATABLE :: sin_grid1_center_lat 549 REAL*8,DIMENSION(:),ALLOCATABLE :: sin_grid1_center_lon 550 INTEGER :: i 551 ! 552 grid2_mask = .TRUE. 553 ! 554 ! 555 nmap = 1 556 ! 557 !*** 558 !*** loop over destination grid 559 !*** 560 ! print*,'grid2_size =',grid2_size 561 ! 562 lastsrc_add=1 563 ! 564 Allocate(cos_grid1_center_lat(size(grid1_center_lat))) 565 Allocate(sin_grid1_center_lat(size(grid1_center_lat))) 566 Allocate(cos_grid1_center_lon(size(grid1_center_lon))) 567 Allocate(sin_grid1_center_lon(size(grid1_center_lon))) 568 569 do i=1,size(grid1_center_lat) 570 cos_grid1_center_lat(i)=cos(grid1_center_lat(i)) 571 sin_grid1_center_lat(i)=sin(grid1_center_lat(i)) 572 ENDDO 573 574 do i=1,size(grid1_center_lon) 575 cos_grid1_center_lon(i)=cos(grid1_center_lon(i)) 576 sin_grid1_center_lon(i)=sin(grid1_center_lon(i)) 577 ENDDO 578 579 grid_loop1: DO dst_add = 1, grid2_size 580 581 IF (.NOT. grid2_mask(dst_add)) CYCLE grid_loop1 582 ! 583 plat = grid2_center_lat(dst_add) 584 plon = grid2_center_lon(dst_add) 585 !*** 586 !*** find nearest square of grid points on source grid 587 !*** 588 CALL grid_search_bilin(src_add, src_lats, src_lons, & 589 plat, plon, grid1_dims, & 590 grid1_center_lat, cos_grid1_center_lat, sin_grid1_center_lat, & 591 grid1_center_lon, cos_grid1_center_lon, sin_grid1_center_lon, & 592 grid1_bound_box, bin_addr1, bin_addr2,lastsrc_add) 593 !*** 594 !*** check to see if points are land points 595 !*** 596 ! 597 IF (src_add(1) > 0) THEN 598 ! 599 DO n=1,4 600 ! if(.not. grid1_mask(src_add(n))) nbmasked = nbmasked + 1 601 IF(.NOT. grid1_mask(src_add(n))) src_add(1) = 0 602 END DO 603 ! 604 ENDIF 605 ! 606 ! 607 608 !*** 609 !*** if point found, find local i,j coordinates for weights 610 !*** 611 IF (src_add(1) > 0) THEN 612 grid2_frac(dst_add) = one 613 !*** 614 !*** iterate to find i,j for bilinear approximation 615 !*** 616 dth1 = src_lats(2) - src_lats(1) 617 dth2 = src_lats(4) - src_lats(1) 618 dth3 = src_lats(3) - src_lats(2) - dth2 619 620 dph1 = src_lons(2) - src_lons(1) 621 dph2 = src_lons(4) - src_lons(1) 622 dph3 = src_lons(3) - src_lons(2) 623 624 IF (dph1 > three*pih) dph1 = dph1 - pi2 625 IF (dph2 > three*pih) dph2 = dph2 - pi2 626 IF (dph3 > three*pih) dph3 = dph3 - pi2 627 IF (dph1 < -three*pih) dph1 = dph1 + pi2 628 IF (dph2 < -three*pih) dph2 = dph2 + pi2 629 IF (dph3 < -three*pih) dph3 = dph3 + pi2 630 631 dph3 = dph3 - dph2 632 633 iguess = half 634 jguess = half 635 636 iter_loop1: DO iter=1,max_iter 637 638 dthp = plat - src_lats(1) - dth1*iguess - & 639 dth2*jguess - dth3*iguess*jguess 640 dphp = plon - src_lons(1) 641 642 IF (dphp > three*pih) dphp = dphp - pi2 643 IF (dphp < -three*pih) dphp = dphp + pi2 644 645 dphp = dphp - dph1*iguess - dph2*jguess - & 646 dph3*iguess*jguess 647 648 mat1 = dth1 + dth3*jguess 649 mat2 = dth2 + dth3*iguess 650 mat3 = dph1 + dph3*jguess 651 mat4 = dph2 + dph3*iguess 652 653 determinant = mat1*mat4 - mat2*mat3 654 655 deli = (dthp*mat4 - mat2*dphp)/determinant 656 delj = (mat1*dphp - dthp*mat3)/determinant 657 658 IF (ABS(deli) < converge .AND. & 659 ABS(delj) < converge) EXIT iter_loop1 660 661 iguess = iguess + deli 662 jguess = jguess + delj 663 664 END DO iter_loop1 665 666 IF (iter <= max_iter) THEN 667 668 !*** 669 !*** successfully found i,j - compute weights 670 !*** 671 672 wgts(1) = (one-iguess)*(one-jguess) 673 wgts(2) = iguess*(one-jguess) 674 wgts(3) = iguess*jguess 675 wgts(4) = (one-iguess)*jguess 676 ! 677 ! 678 CALL store_link_bilin(dst_add, src_add, wgts, nmap) 679 680 ELSE 681 PRINT *,'Point coords: ',plat,plon 682 PRINT *,'Dest grid lats: ',src_lats 683 PRINT *,'Dest grid lons: ',src_lons 684 PRINT *,'Dest grid addresses: ',src_add 685 PRINT *,'Current i,j : ',iguess, jguess 686 STOP 'Iteration for i,j exceed max iteration count' 687 ENDIF 688 ! 689 !*** 690 !*** search for bilinear failed - use a distance-weighted 691 !*** average instead (this is typically near the pole) 692 !*** 693 ELSE IF (src_add(1) < 0) THEN 694 695 src_add = ABS(src_add) 696 icount = 0 697 DO n=1,4 698 ! 699 IF (grid1_mask(src_add(n))) THEN 700 icount = icount + 1 701 ELSE 702 src_lats(n) = zero 703 ENDIF 704 ! 705 END DO 706 707 IF (icount > 0) THEN 708 ! 709 !*** renormalize weights 710 ! 711 sum_wgts = SUM(src_lats) 712 wgts(1) = src_lats(1)/sum_wgts 713 wgts(2) = src_lats(2)/sum_wgts 714 wgts(3) = src_lats(3)/sum_wgts 715 wgts(4) = src_lats(4)/sum_wgts 716 ! 717 grid2_frac(dst_add) = one 718 CALL store_link_bilin(dst_add, src_add, wgts, nmap) 719 ENDIF 720 721 ! 722 ENDIF 723 END DO grid_loop1 724 ! 725 ! Call sort_add(grid2_add_map1, grid1_add_map1, wts_map1) 726 ! 727 ! 728 !----------------------------------------------------------------------- 729 730 deallocate(cos_grid1_center_lat) 731 deallocate(cos_grid1_center_lon) 732 deallocate(sin_grid1_center_lat) 733 deallocate(sin_grid1_center_lon) 734 END SUBROUTINE remap_bilin 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 !*********************************************************************** 754 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 755 !************************************************************************ 756 ! SUBROUTINE GRID_SEARCH_BILIN 757 !************************************************************************ 758 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 759 ! 760 SUBROUTINE grid_search_bilin(src_add, src_lats, src_lons, & 761 plat, plon, src_grid_dims, & 762 src_center_lat, cos_src_center_lat, sin_src_center_lat, & 763 src_center_lon, cos_src_center_lon, sin_src_center_lon,& 764 src_grid_bound_box, & 765 src_bin_add, dst_bin_add,lastsrc_add) 766 767 !----------------------------------------------------------------------- 768 ! 769 ! this routine finds the location of the search point plat, plon 770 ! in the source grid and returns the corners needed for a bilinear 771 ! interpolation. 772 ! 773 !----------------------------------------------------------------------- 774 775 !----------------------------------------------------------------------- 776 ! output variables 777 !----------------------------------------------------------------------- 778 ! 779 ! address of each corner point enclosing P 780 ! 781 INTEGER,DIMENSION(4) :: src_add 782 REAL*8,DIMENSION(4) :: src_lats,src_lons 783 ! 784 !----------------------------------------------------------------------- 785 ! input variables 786 !----------------------------------------------------------------------- 787 ! 788 ! latitude, longitude of the search point 789 ! 790 REAL*8, INTENT(in) :: plat,plon 791 ! 792 ! size of each src grid dimension 793 ! 794 INTEGER, DIMENSION(2), INTENT(in) :: src_grid_dims 795 ! 796 ! latitude, longitude of each src grid center 797 ! 798 REAL*8, DIMENSION(:), INTENT(in) :: src_center_lat,src_center_lon 799 REAL*8, DIMENSION(:), INTENT(in) :: cos_src_center_lat 800 REAL*8, DIMENSION(:), INTENT(in) :: sin_src_center_lat 801 REAL*8, DIMENSION(:), INTENT(in) :: cos_src_center_lon 802 REAL*8, DIMENSION(:), INTENT(in) :: sin_src_center_lon 803 ! 804 ! bound box for source grid 805 ! 806 REAL*8, DIMENSION(:,:), INTENT(in) :: src_grid_bound_box 807 ! 808 ! latitude bins for restricting searches 809 ! 810 INTEGER, DIMENSION(:,:), INTENT(in) ::src_bin_add,dst_bin_add 811 812 INTEGER,OPTIONAL :: lastsrc_add 813 INTEGER :: loopsrc,l1,l2 814 ! 815 !----------------------------------------------------------------------- 816 ! local variables 817 !----------------------------------------------------------------------- 818 ! 819 INTEGER :: n,next_n,srch_add,nx, ny,min_add, max_add, & 820 i, j, jp1, ip1, n_add, e_add, ne_add 821 822 823 REAL*8 :: vec1_lat, vec1_lon,vec2_lat, vec2_lon, cross_product, & 824 cross_product_last,coslat_dst, sinlat_dst, coslon_dst, & 825 sinlon_dst,dist_min, distance 826 827 !----------------------------------------------------------------------- 828 ! restrict search first using bins 829 !----------------------------------------------------------------------- 830 831 src_add = 0 832 833 min_add = SIZE(src_center_lat) 834 max_add = 1 835 DO n=1,num_srch_bins 836 IF (plat >= bin_lats(1,n) .AND. plat <= bin_lats(2,n) .AND. & 837 plon >= bin_lons(1,n) .AND. plon <= bin_lons(2,n)) THEN 838 min_add = MIN(min_add, src_bin_add(1,n)) 839 max_add = MAX(max_add, src_bin_add(2,n)) 840 ENDIF 841 END DO 842 843 !----------------------------------------------------------------------- 844 ! now perform a more detailed search 845 !----------------------------------------------------------------------- 846 847 nx = src_grid_dims(1) 848 ny = src_grid_dims(2) 849 850 loopsrc=0 851 DO WHILE (loopsrc <= max_add) 852 853 854 l1=MAX(min_add,lastsrc_add-loopsrc) 855 l2=MIN(max_add,lastsrc_add+loopsrc) 856 857 loopsrc = loopsrc+1 858 859 srch_loop: DO srch_add = l1,l2,MAX(l2-l1,1) 860 861 !*** first check bounding box 862 863 IF (plat <= src_grid_bound_box(2,srch_add) .AND. & 864 plat >= src_grid_bound_box(1,srch_add) .AND. & 865 plon <= src_grid_bound_box(4,srch_add) .AND. & 866 plon >= src_grid_bound_box(3,srch_add)) THEN 867 !*** 868 !*** we are within bounding box so get really serious 869 !*** 870 !*** determine neighbor addresses 871 ! 872 j = (srch_add - 1)/nx +1 873 i = srch_add - (j-1)*nx 874 ! 875 IF (i < nx) THEN 876 ip1 = i + 1 877 ELSE 878 ip1 = 1 879 ENDIF 880 ! 881 IF (j < ny) THEN 882 jp1 = j+1 883 ELSE 884 jp1 = 1 885 ENDIF 886 ! 887 n_add = (jp1 - 1)*nx + i 888 e_add = (j - 1)*nx + ip1 889 ne_add = (jp1 - 1)*nx + ip1 890 ! 891 src_lats(1) = src_center_lat(srch_add) 892 src_lats(2) = src_center_lat(e_add) 893 src_lats(3) = src_center_lat(ne_add) 894 src_lats(4) = src_center_lat(n_add) 895 ! 896 src_lons(1) = src_center_lon(srch_add) 897 src_lons(2) = src_center_lon(e_add) 898 src_lons(3) = src_center_lon(ne_add) 899 src_lons(4) = src_center_lon(n_add) 900 ! 901 !*** 902 !*** for consistency, we must make sure all lons are in 903 !*** same 2pi interval 904 !*** 905 ! 906 vec1_lon = src_lons(1) - plon 907 IF (vec1_lon > pi) THEN 908 src_lons(1) = src_lons(1) - pi2 909 ELSE IF (vec1_lon < -pi) THEN 910 src_lons(1) = src_lons(1) + pi2 911 ENDIF 912 DO n=2,4 913 vec1_lon = src_lons(n) - src_lons(1) 914 IF (vec1_lon > pi) THEN 915 src_lons(n) = src_lons(n) - pi2 916 ELSE IF (vec1_lon < -pi) THEN 917 src_lons(n) = src_lons(n) + pi2 918 ENDIF 919 END DO 920 ! 921 corner_loop: DO n=1,4 922 next_n = MOD(n,4) + 1 923 !*** 924 !*** here we take the cross product of the vector making 925 !*** up each box side with the vector formed by the vertex 926 !*** and search point. if all the cross products are 927 !*** positive, the point is contained in the box. 928 !*** 929 vec1_lat = src_lats(next_n) - src_lats(n) 930 vec1_lon = src_lons(next_n) - src_lons(n) 931 vec2_lat = plat - src_lats(n) 932 vec2_lon = plon - src_lons(n) 933 !*** 934 !*** check for 0,2pi crossings 935 !*** 936 IF (vec1_lon > three*pih) THEN 937 vec1_lon = vec1_lon - pi2 938 ELSE IF (vec1_lon < -three*pih) THEN 939 vec1_lon = vec1_lon + pi2 940 ENDIF 941 IF (vec2_lon > three*pih) THEN 942 vec2_lon = vec2_lon - pi2 943 ELSE IF (vec2_lon < -three*pih) THEN 944 vec2_lon = vec2_lon + pi2 945 ENDIF 946 ! 947 cross_product = vec1_lon*vec2_lat - vec2_lon*vec1_lat 948 ! 949 !*** 950 !*** if cross product is less than zero, this cell 951 !*** doesn't work 952 !*** 953 IF (n == 1) cross_product_last = cross_product 954 IF (cross_product*cross_product_last < zero) & 955 EXIT corner_loop 956 cross_product_last = cross_product 957 ! 958 END DO corner_loop 959 !*** 960 !*** if cross products all same sign, we found the location 961 !*** 962 IF (n > 4) THEN 963 src_add(1) = srch_add 964 src_add(2) = e_add 965 src_add(3) = ne_add 966 src_add(4) = n_add 967 ! 968 lastsrc_add = srch_add 969 RETURN 970 ENDIF 971 !*** 972 !*** otherwise move on to next cell 973 !*** 974 ENDIF !bounding box check 975 END DO srch_loop 976 977 978 ENDDO 979 980 981 !*** 982 !*** if no cell found, point is likely either in a box that 983 !*** straddles either pole or is outside the grid. fall back 984 !*** to a distance-weighted average of the four closest 985 !*** points. go ahead and compute weights here, but store 986 !*** in src_lats and return -add to prevent the parent 987 !*** routine from computing bilinear weights 988 !*** 989 !print *,'Could not find location for ',plat,plon 990 !print *,'Using nearest-neighbor average for this point' 991 ! 992 coslat_dst = COS(plat) 993 sinlat_dst = SIN(plat) 994 coslon_dst = COS(plon) 995 sinlon_dst = SIN(plon) 996 ! 997 dist_min = bignum 998 src_lats = bignum 999 DO srch_add = min_add,max_add 1000 ! distance = ACOS(coslat_dst*COS(src_center_lat(srch_add))* & 1001 ! (coslon_dst*COS(src_center_lon(srch_add)) + & 1002 ! sinlon_dst*SIN(src_center_lon(srch_add)))+ & 1003 ! sinlat_dst*SIN(src_center_lat(srch_add))) 1004 1005 distance = ACOS(coslat_dst*cos_src_center_lat(srch_add)* & 1006 (coslon_dst*cos_src_center_lon(srch_add) + & 1007 sinlon_dst*sin_src_center_lon(srch_add))+ & 1008 sinlat_dst*sin_src_center_lat(srch_add)) 1009 1010 IF (distance < dist_min) THEN 1011 sort_loop: DO n=1,4 1012 IF (distance < src_lats(n)) THEN 1013 DO i=4,n+1,-1 1014 src_add (i) = src_add (i-1) 1015 src_lats(i) = src_lats(i-1) 1016 END DO 1017 src_add (n) = -srch_add 1018 src_lats(n) = distance 1019 dist_min = src_lats(4) 1020 EXIT sort_loop 1021 ENDIF 1022 END DO sort_loop 1023 ENDIF 1024 END DO 1025 ! 1026 src_lons = one/(src_lats + tiny) 1027 distance = SUM(src_lons) 1028 src_lats = src_lons/distance 1029 1030 !----------------------------------------------------------------------- 1031 1032 END SUBROUTINE grid_search_bilin 1033 1034 1035 !*********************************************************************** 1036 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1037 !************************************************************************ 1038 ! SUBROUTINE STORE_LINK_BILIN 1039 !************************************************************************ 1040 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1041 1042 SUBROUTINE store_link_bilin(dst_add, src_add, weights, nmap) 1043 1044 !----------------------------------------------------------------------- 1045 ! this routine stores the address and weight for four links 1046 ! associated with one destination point in the appropriate address 1047 ! and weight arrays and resizes those arrays if necessary. 1048 !----------------------------------------------------------------------- 1049 1050 !----------------------------------------------------------------------- 1051 ! input variables 1052 !----------------------------------------------------------------------- 1053 ! 1054 INTEGER :: dst_add,nmap 1055 ! 1056 INTEGER, DIMENSION(4) :: src_add 1057 ! 1058 REAL*8, DIMENSION(4) :: weights 1059 1060 !----------------------------------------------------------------------- 1061 ! 1062 ! local variables 1063 ! 1064 !----------------------------------------------------------------------- 1065 1066 INTEGER :: n,num_links_old 1067 1068 !----------------------------------------------------------------------- 1069 ! increment number of links and check to see if remap arrays need 1070 ! to be increased to accomodate the new link. then store the 1071 ! link. 1072 !----------------------------------------------------------------------- 1073 1074 num_links_old = num_links_map1 1075 num_links_map1 = num_links_old + 4 1076 1077 IF (num_links_map1 > max_links_map1) & 1078 CALL resize_remap_vars(1,resize_increment) 1079 1080 DO n=1,4 1081 grid1_add_map1(num_links_old+n) = src_add(n) 1082 grid2_add_map1(num_links_old+n) = dst_add 1083 wts_map1 (1,num_links_old+n) = weights(n) 1084 END DO 1085 1086 !----------------------------------------------------------------------- 1087 1088 END SUBROUTINE store_link_bilin 1089 1090 1091 1092 1093 1094 1095 1096 1097 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1098 !************************************************************************ 1099 ! SUBROUTINE INIT_REMAP_VARS 1100 !************************************************************************ 1101 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1102 ! 1103 SUBROUTINE init_remap_vars 1104 1105 !----------------------------------------------------------------------- 1106 ! 1107 ! this routine initializes some variables and provides an initial 1108 ! allocation of arrays (fairly large so frequent resizing 1109 ! unnecessary). 1110 ! 1111 !----------------------------------------------------------------------- 1112 1113 !----------------------------------------------------------------------- 1114 ! determine the number of weights 1115 !----------------------------------------------------------------------- 1116 num_wts = 1 ! bilinear interpolation 1117 !----------------------------------------------------------------------- 1118 ! initialize num_links and set max_links to four times the largest 1119 ! of the destination grid sizes initially (can be changed later). 1120 ! set a default resize increment to increase the size of link 1121 ! arrays if the number of links exceeds the initial size 1122 !----------------------------------------------------------------------- 1123 1124 num_links_map1 = 0 1125 max_links_map1 = 4*grid2_size 1126 IF (num_maps > 1) THEN 1127 num_links_map2 = 0 1128 max_links_map1 = MAX(4*grid1_size,4*grid2_size) 1129 max_links_map2 = max_links_map1 1130 ENDIF 1131 1132 resize_increment = 0.1*MAX(grid1_size,grid2_size) 1133 1134 !----------------------------------------------------------------------- 1135 ! allocate address and weight arrays for mapping 1 1136 !----------------------------------------------------------------------- 1137 1138 ALLOCATE (grid1_add_map1(max_links_map1), & 1139 grid2_add_map1(max_links_map1), & 1140 wts_map1(num_wts, max_links_map1)) 1141 1142 grid1_add_map1 = 0. 1143 grid2_add_map1 = 0. 1144 wts_map1 = 0. 1145 1146 !----------------------------------------------------------------------- 1147 1148 END SUBROUTINE init_remap_vars 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 !*********************************************************************** 1164 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1165 !************************************************************************ 1166 ! SUBROUTINE RESIZE_REMAP_VAR 1167 !************************************************************************ 1168 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1169 1170 SUBROUTINE resize_remap_vars(nmap, increment) 1171 1172 !----------------------------------------------------------------------- 1173 ! this routine resizes remapping arrays by increasing(decreasing) 1174 ! the max_links by increment 1175 !----------------------------------------------------------------------- 1176 1177 !----------------------------------------------------------------------- 1178 ! input variables 1179 !----------------------------------------------------------------------- 1180 1181 INTEGER :: & 1182 nmap, & ! identifies which mapping array to resize 1183 increment ! the number of links to add(subtract) to arrays 1184 1185 !----------------------------------------------------------------------- 1186 ! local variables 1187 !----------------------------------------------------------------------- 1188 1189 INTEGER :: & 1190 ierr, & ! error flag 1191 mxlinks ! size of link arrays 1192 1193 INTEGER, DIMENSION(:), ALLOCATABLE :: & 1194 add1_tmp, & ! temp array for resizing address arrays 1195 add2_tmp ! temp array for resizing address arrays 1196 ! 1197 ! temp array for resizing weight arrays 1198 ! 1199 REAL*8, DIMENSION(:,:), ALLOCATABLE :: wts_tmp 1200 ! 1201 !----------------------------------------------------------------------- 1202 !*** 1203 !*** allocate temporaries to hold original values 1204 !*** 1205 mxlinks = SIZE(grid1_add_map1) 1206 ALLOCATE (add1_tmp(mxlinks), add2_tmp(mxlinks), & 1207 wts_tmp(num_wts,mxlinks)) 1208 1209 add1_tmp = grid1_add_map1 1210 add2_tmp = grid2_add_map1 1211 wts_tmp = wts_map1 1212 1213 !*** 1214 !*** deallocate originals and increment max_links then 1215 !*** reallocate arrays at new size 1216 !*** 1217 1218 DEALLOCATE (grid1_add_map1, grid2_add_map1, wts_map1) 1219 max_links_map1 = mxlinks + increment 1220 ALLOCATE (grid1_add_map1(max_links_map1), & 1221 grid2_add_map1(max_links_map1), & 1222 wts_map1(num_wts,max_links_map1)) 1223 !*** 1224 !*** restore original values from temp arrays and 1225 !*** deallocate temps 1226 !*** 1227 mxlinks = MIN(mxlinks, max_links_map1) 1228 grid1_add_map1(1:mxlinks) = add1_tmp (1:mxlinks) 1229 grid2_add_map1(1:mxlinks) = add2_tmp (1:mxlinks) 1230 wts_map1 (:,1:mxlinks) = wts_tmp(:,1:mxlinks) 1231 DEALLOCATE(add1_tmp, add2_tmp, wts_tmp) 1232 1233 !----------------------------------------------------------------------- 1234 ! 1235 END SUBROUTINE resize_remap_vars 1236 ! 1237 !************************************************************************ 1238 ! 1239 END MODULE bilinear_interp 1240 1241 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1242
Note: See TracChangeset
for help on using the changeset viewer.
