source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/TOOLS/WEIGHTS/src/remap_bicubic.f90 @ 5443

Last change on this file since 5443 was 5443, checked in by davestorkey, 5 years ago

Update 2015/dev_r5021_UKMO1_CICE_coupling branch to revision 5442 of the trunk.

File size: 28.0 KB
Line 
1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2!
3!     this module contains necessary routines for performing an
4!     bicubic interpolation.
5!
6!-----------------------------------------------------------------------
7!
8!     CVS:$Id$
9!
10!     Copyright (c) 1997, 1998 the Regents of the University of
11!       California.
12!
13!     This software and ancillary information (herein called software)
14!     called SCRIP is made available under the terms described here. 
15!     The software has been approved for release with associated
16!     LA-CC Number 98-45.
17!
18!     Unless otherwise indicated, this software has been authored
19!     by an employee or employees of the University of California,
20!     operator of the Los Alamos National Laboratory under Contract
21!     No. W-7405-ENG-36 with the U.S. Department of Energy.  The U.S.
22!     Government has rights to use, reproduce, and distribute this
23!     software.  The public may copy and use this software without
24!     charge, provided that this Notice and any statement of authorship
25!     are reproduced on all copies.  Neither the Government nor the
26!     University makes any warranty, express or implied, or assumes
27!     any liability or responsibility for the use of this software.
28!
29!     If software is modified to produce derivative works, such modified
30!     software should be clearly marked, so as not to confuse it with
31!     the version available from Los Alamos National Laboratory.
32!
33!***********************************************************************
34
35      module remap_bicubic
36
37!-----------------------------------------------------------------------
38
39      use kinds_mod     ! defines common data types
40      use constants     ! defines common constants
41      use grids         ! module containing grid info
42      use remap_vars    ! module containing remap info
43
44      implicit none
45
46!-----------------------------------------------------------------------
47
48      integer (kind=int_kind), parameter :: &
49          max_iter = 100   ! max iteration count for i,j iteration
50
51      real (kind=dbl_kind), parameter :: &
52           converge = 1.e-10_dbl_kind  ! convergence criterion
53
54!***********************************************************************
55
56      contains
57
58!***********************************************************************
59
60      subroutine remap_bicub
61
62!-----------------------------------------------------------------------
63!
64!     this routine computes the weights for a bicubic interpolation.
65!
66!-----------------------------------------------------------------------
67
68!-----------------------------------------------------------------------
69!
70!     local variables
71!
72!-----------------------------------------------------------------------
73
74      integer (kind=int_kind) :: n,icount, &
75           dst_add,         & ! destination address
76           iter,            & ! iteration counter
77           nmap            ! index of current map being computed
78
79      integer (kind=int_kind), dimension(4) ::  &
80           src_add         ! address for the four source points
81
82      real (kind=dbl_kind), dimension(4)  :: &
83           src_lats,        & ! latitudes  of four bilinear corners
84           src_lons        ! longitudes of four bilinear corners
85
86      real (kind=dbl_kind), dimension(4,4)  :: &
87           wgts            ! bicubic weights for four corners
88
89      real (kind=dbl_kind) :: &
90           plat, plon,        & ! lat/lon coords of destination point
91           iguess, jguess,    & ! current guess for bilinear coordinate
92           thguess, phguess,  & ! current guess for lat/lon coordinate
93           deli, delj,        & ! corrections to i,j
94           dth1, dth2, dth3,  & ! some latitude  differences
95           dph1, dph2, dph3,  & ! some longitude differences
96           dthp, dphp,        & ! difference between point and sw corner
97           mat1, mat2, mat3, mat4,  & ! matrix elements
98           determinant,       & ! matrix determinant
99           sum_wgts,          & ! sum of weights for normalization
100           w1,w2,w3,w4,w5,w6,w7,w8,  & ! 16 bicubic weight functions
101           w9,w10,w11,w12,w13,w14,w15,w16
102
103!-----------------------------------------------------------------------
104!
105!     compute mappings from grid1 to grid2
106!
107!-----------------------------------------------------------------------
108
109      nmap = 1
110      if (grid1_rank /= 2) then
111        stop 'Can not do bicubic interpolation when grid_rank /= 2'
112      endif
113
114      !***
115      !*** loop over destination grid
116      !***
117
118      grid_loop1: do dst_add = 1, grid2_size
119
120        if (.not. grid2_mask(dst_add)) cycle grid_loop1
121
122        plat = grid2_center_lat(dst_add)
123        plon = grid2_center_lon(dst_add)
124
125!-----------------------------------------------------------------------
126!
127!       find nearest square of grid points on source grid
128!
129!-----------------------------------------------------------------------
130
131        call grid_search_bicub(src_add, src_lats, src_lons,  &
132                               plat, plon, grid1_dims, &
133                               grid1_center_lat, grid1_center_lon, &
134                               grid1_bound_box, bin_addr1, bin_addr2)
135
136        !***
137        !*** check to see if points are land points
138        !***
139
140        if (src_add(1) > 0) then
141          do n=1,4
142            if (.not. grid1_mask(src_add(n))) src_add(1) = 0
143          end do
144        endif
145
146!-----------------------------------------------------------------------
147!
148!       if point found, find local i,j coordinates for weights
149!
150!-----------------------------------------------------------------------
151
152        if (src_add(1) > 0) then
153
154          grid2_frac(dst_add) = one
155
156          !***
157          !*** iterate to find i,j for bicubic approximation
158          !***
159
160          dth1 = src_lats(2) - src_lats(1)
161          dth2 = src_lats(4) - src_lats(1)
162          dth3 = src_lats(3) - src_lats(2) - dth2
163
164          dph1 = src_lons(2) - src_lons(1)
165          dph2 = src_lons(4) - src_lons(1)
166          dph3 = src_lons(3) - src_lons(2)
167
168          if (dph1 >  three*pih) dph1 = dph1 - pi2
169          if (dph2 >  three*pih) dph2 = dph2 - pi2
170          if (dph3 >  three*pih) dph3 = dph3 - pi2
171          if (dph1 < -three*pih) dph1 = dph1 + pi2
172          if (dph2 < -three*pih) dph2 = dph2 + pi2
173          if (dph3 < -three*pih) dph3 = dph3 + pi2
174
175          dph3 = dph3 - dph2
176
177          iguess = half
178          jguess = half
179
180          iter_loop1: do iter=1,max_iter
181
182            dthp = plat - src_lats(1) - dth1*iguess - &
183                          dth2*jguess - dth3*iguess*jguess
184            dphp = plon - src_lons(1)
185
186            if (dphp >  three*pih) dphp = dphp - pi2
187            if (dphp < -three*pih) dphp = dphp + pi2
188
189            dphp = dphp - dph1*iguess - dph2*jguess -  &
190                          dph3*iguess*jguess
191
192            mat1 = dth1 + dth3*jguess
193            mat2 = dth2 + dth3*iguess
194            mat3 = dph1 + dph3*jguess
195            mat4 = dph2 + dph3*iguess
196
197            determinant = mat1*mat4 - mat2*mat3
198
199            deli = (dthp*mat4 - mat2*dphp)/determinant
200            delj = (mat1*dphp - dthp*mat3)/determinant
201
202            if (abs(deli) < converge .and.  &
203                abs(delj) < converge) exit iter_loop1
204
205            iguess = iguess + deli
206            jguess = jguess + delj
207
208          end do iter_loop1
209
210          if (iter <= max_iter) then
211
212!-----------------------------------------------------------------------
213!
214!           successfully found i,j - compute weights
215!
216!-----------------------------------------------------------------------
217
218            wgts(1,1) = (one - jguess**2*(three-two*jguess))* &
219                        (one - iguess**2*(three-two*iguess))
220            wgts(1,2) = (one - jguess**2*(three-two*jguess))* &
221                               iguess**2*(three-two*iguess)
222            wgts(1,3) =        jguess**2*(three-two*jguess)* &
223                               iguess**2*(three-two*iguess)
224            wgts(1,4) =        jguess**2*(three-two*jguess)* &
225                        (one - iguess**2*(three-two*iguess))
226            wgts(2,1) = (one - jguess**2*(three-two*jguess))* &
227                               iguess*(iguess-one)**2
228            wgts(2,2) = (one - jguess**2*(three-two*jguess))* &
229                               iguess**2*(iguess-one)
230            wgts(2,3) =        jguess**2*(three-two*jguess)* &
231                               iguess**2*(iguess-one)
232            wgts(2,4) =        jguess**2*(three-two*jguess)* &
233                               iguess*(iguess-one)**2
234            wgts(3,1) =        jguess*(jguess-one)**2* &
235                        (one - iguess**2*(three-two*iguess))
236            wgts(3,2) =        jguess*(jguess-one)**2* &
237                               iguess**2*(three-two*iguess)
238            wgts(3,3) =        jguess**2*(jguess-one)* &
239                               iguess**2*(three-two*iguess)
240            wgts(3,4) =        jguess**2*(jguess-one)* &
241                        (one - iguess**2*(three-two*iguess))
242            wgts(4,1) =        iguess*(iguess-one)**2* &
243                               jguess*(jguess-one)**2
244            wgts(4,2) =        iguess**2*(iguess-one)* &
245                               jguess*(jguess-one)**2
246            wgts(4,3) =        iguess**2*(iguess-one)* &
247                               jguess**2*(jguess-one)
248            wgts(4,4) =        iguess*(iguess-one)**2* &
249                               jguess**2*(jguess-one)
250
251            call store_link_bicub(dst_add, src_add, wgts, nmap)
252
253          else
254            stop 'Iteration for i,j exceed max iteration count'
255          endif
256
257!-----------------------------------------------------------------------
258!
259!       search for bilinear failed - use a distance-weighted
260!       average instead (this is typically near the pole)
261!
262!-----------------------------------------------------------------------
263
264        else if (src_add(1) < 0) then
265
266          src_add = abs(src_add)
267
268          icount = 0
269          do n=1,4
270            if (grid1_mask(src_add(n))) then
271              icount = icount + 1
272            else
273              src_lats(n) = zero
274            endif
275          end do
276
277          if (icount > 0) then
278            !*** renormalize weights
279
280            sum_wgts = sum(src_lats)
281            wgts(1,1) = src_lats(1)/sum_wgts
282            wgts(1,2) = src_lats(2)/sum_wgts
283            wgts(1,3) = src_lats(3)/sum_wgts
284            wgts(1,4) = src_lats(4)/sum_wgts
285            wgts(2:4,:) = zero
286
287            grid2_frac(dst_add) = one
288            call store_link_bicub(dst_add, src_add, wgts, nmap)
289          endif
290
291        endif
292      end do grid_loop1
293
294!-----------------------------------------------------------------------
295!
296!     compute mappings from grid2 to grid1 if necessary
297!
298!-----------------------------------------------------------------------
299
300      if (num_maps > 1) then
301
302      nmap = 2
303      if (grid2_rank /= 2) then
304        stop 'Can not do bicubic interpolation when grid_rank /= 2'
305      endif
306
307      !***
308      !*** loop over destination grid
309      !***
310
311      grid_loop2: do dst_add = 1, grid1_size
312
313        if (.not. grid1_mask(dst_add)) cycle grid_loop2
314
315        plat = grid1_center_lat(dst_add)
316        plon = grid1_center_lon(dst_add)
317
318        !***
319        !*** find nearest square of grid points on source grid
320        !***
321
322        call grid_search_bicub(src_add, src_lats, src_lons,  &
323                               plat, plon, grid2_dims, &
324                               grid2_center_lat, grid2_center_lon, &
325                               grid2_bound_box, bin_addr2, bin_addr1)
326
327        !***
328        !*** check to see if points are land points
329        !***
330
331        if (src_add(1) > 0) then
332          do n=1,4
333            if (.not. grid2_mask(src_add(n))) src_add(1) = 0
334          end do
335        endif
336
337        !***
338        !*** if point found, find i,j coordinates for weights
339        !***
340
341        if (src_add(1) > 0) then
342
343          grid1_frac(dst_add) = one
344
345          !***
346          !*** iterate to find i,j for bilinear approximation
347          !***
348
349          dth1 = src_lats(2) - src_lats(1)
350          dth2 = src_lats(4) - src_lats(1)
351          dth3 = src_lats(3) - src_lats(2) - dth2
352
353          dph1 = src_lons(2) - src_lons(1)
354          dph2 = src_lons(4) - src_lons(1)
355          dph3 = src_lons(3) - src_lons(2)
356
357          if (dph1 >  pi) dph1 = dph1 - pi2
358          if (dph2 >  pi) dph2 = dph2 - pi2
359          if (dph3 >  pi) dph3 = dph3 - pi2
360          if (dph1 < -pi) dph1 = dph1 + pi2
361          if (dph2 < -pi) dph2 = dph2 + pi2
362          if (dph3 < -pi) dph3 = dph3 + pi2
363
364          dph3 = dph3 - dph2
365
366          iguess = zero
367          jguess = zero
368
369          iter_loop2: do iter=1,max_iter
370
371            dthp = plat - src_lats(1) - dth1*iguess - &
372                          dth2*jguess - dth3*iguess*jguess
373            dphp = plon - src_lons(1)
374
375            if (dphp >  pi) dphp = dphp - pi2
376            if (dphp < -pi) dphp = dphp + pi2
377
378            dphp = dphp - dph1*iguess - dph2*jguess -  &
379                          dph3*iguess*jguess
380
381            mat1 = dth1 + dth3*jguess
382            mat2 = dth2 + dth3*iguess
383            mat3 = dph1 + dph3*jguess
384            mat4 = dph2 + dph3*iguess
385
386            determinant = mat1*mat4 - mat2*mat3
387
388            deli = (dthp*mat4 - mat2*dphp)/determinant
389            delj = (mat1*dphp - dthp*mat3)/determinant
390
391            if (abs(deli) < converge .and.  &
392                abs(delj) < converge) exit iter_loop2
393
394            iguess = iguess + deli
395            jguess = jguess + delj
396
397          end do iter_loop2
398
399          if (iter <= max_iter) then
400
401            !***
402            !*** successfully found i,j - compute weights
403            !***
404
405            wgts(1,1) = (one - jguess**2*(three-two*jguess))* &
406                        (one - iguess**2*(three-two*iguess))
407            wgts(1,2) = (one - jguess**2*(three-two*jguess))* &
408                               iguess**2*(three-two*iguess)
409            wgts(1,3) =        jguess**2*(three-two*jguess)* &
410                               iguess**2*(three-two*iguess)
411            wgts(1,4) =        jguess**2*(three-two*jguess)* &
412                        (one - iguess**2*(three-two*iguess))
413            wgts(2,1) = (one - jguess**2*(three-two*jguess))* &
414                               iguess*(iguess-one)**2
415            wgts(2,2) = (one - jguess**2*(three-two*jguess))* &
416                               iguess**2*(iguess-one)
417            wgts(2,3) =        jguess**2*(three-two*jguess)* &
418                               iguess**2*(iguess-one)
419            wgts(2,4) =        jguess**2*(three-two*jguess)* &
420                               iguess*(iguess-one)**2
421            wgts(3,1) =        jguess*(jguess-one)**2* &
422                        (one - iguess**2*(three-two*iguess))
423            wgts(3,2) =        jguess*(jguess-one)**2* &
424                               iguess**2*(three-two*iguess)
425            wgts(3,3) =        jguess**2*(jguess-one)* &
426                               iguess**2*(three-two*iguess)
427            wgts(3,4) =        jguess**2*(jguess-one)* &
428                        (one - iguess**2*(three-two*iguess))
429            wgts(4,1) =        iguess*(iguess-one)**2* &
430                               jguess*(jguess-one)**2
431            wgts(4,2) =        iguess**2*(iguess-one)* &
432                               jguess*(jguess-one)**2
433            wgts(4,3) =        iguess**2*(iguess-one)* &
434                               jguess**2*(jguess-one)
435            wgts(4,4) =        iguess*(iguess-one)**2* &
436                               jguess**2*(jguess-one)
437
438            call store_link_bicub(dst_add, src_add, wgts, nmap)
439
440          else
441            stop 'Iteration for i,j exceed max iteration count'
442          endif
443
444        !***
445        !*** search for bilinear failed - us a distance-weighted
446        !*** average instead
447        !***
448
449        else if (src_add(1) < 0) then
450
451          src_add = abs(src_add)
452
453          icount = 0
454          do n=1,4
455            if (grid2_mask(src_add(n))) then
456              icount = icount + 1
457            else
458              src_lats(n) = zero
459            endif
460          end do
461
462          if (icount > 0) then
463            !*** renormalize weights
464
465            sum_wgts = sum(src_lats)
466            wgts(1,1) = src_lats(1)/sum_wgts
467            wgts(1,2) = src_lats(2)/sum_wgts
468            wgts(1,3) = src_lats(3)/sum_wgts
469            wgts(1,4) = src_lats(4)/sum_wgts
470            wgts(2:4,:) = zero
471
472            grid1_frac(dst_add) = one
473            call store_link_bicub(dst_add, src_add, wgts, nmap)
474          endif
475
476        endif
477      end do grid_loop2
478
479      endif ! nmap=2
480
481!-----------------------------------------------------------------------
482
483      end subroutine remap_bicub
484
485!***********************************************************************
486
487      subroutine grid_search_bicub(src_add, src_lats, src_lons,  &
488                                   plat, plon, src_grid_dims, &
489                                   src_center_lat, src_center_lon, &
490                                   src_bound_box, &
491                                   src_bin_add, dst_bin_add)
492
493!-----------------------------------------------------------------------
494!
495!     this routine finds the location of the search point plat, plon
496!     in the source grid and returns the corners needed for a bicubic
497!     interpolation.
498!
499!-----------------------------------------------------------------------
500
501!-----------------------------------------------------------------------
502!
503!     output variables
504!
505!-----------------------------------------------------------------------
506
507      integer (kind=int_kind), dimension(4), intent(out) :: &
508              src_add  ! address of each corner point enclosing P
509
510      real (kind=dbl_kind), dimension(4), intent(out) :: &
511              src_lats,  & ! latitudes  of the four corner points
512              src_lons  ! longitudes of the four corner points
513
514!-----------------------------------------------------------------------
515!
516!     input variables
517!
518!-----------------------------------------------------------------------
519
520      real (kind=dbl_kind), intent(in) :: &
521              plat,    & ! latitude  of the search point
522              plon    ! longitude of the search point
523
524      integer (kind=int_kind), dimension(2), intent(in) :: &
525              src_grid_dims  ! size of each src grid dimension
526
527      real (kind=dbl_kind), dimension(:), intent(in) :: &
528              src_center_lat,  & ! latitude  of each src grid center
529              src_center_lon  ! longitude of each src grid center
530
531      real (kind=dbl_kind), dimension(:,:), intent(in) :: &
532              src_bound_box   ! bounding box for src grid search
533
534      integer (kind=int_kind), dimension(:,:), intent(in) :: &
535              src_bin_add,     & ! search bins for restricting
536              dst_bin_add     ! searches
537
538!-----------------------------------------------------------------------
539!
540!     local variables
541!
542!-----------------------------------------------------------------------
543
544      integer (kind=int_kind) :: n, next_n, srch_add,    & ! dummy indices
545          nx, ny,             & ! dimensions of src grid
546          min_add, max_add,   & ! addresses for restricting search
547          i, j, jp1, ip1, n_add, e_add, ne_add  ! addresses
548
549      real (kind=dbl_kind) ::   & ! vectors for cross-product check
550            vec1_lat, vec1_lon, &
551            vec2_lat, vec2_lon, cross_product, cross_product_last, &
552            coslat_dst, sinlat_dst, coslon_dst, sinlon_dst, &
553            dist_min, distance ! for computing dist-weighted avg
554
555!-----------------------------------------------------------------------
556!
557!     restrict search first using search bins.
558!
559!-----------------------------------------------------------------------
560
561      src_add = 0
562
563      min_add = size(src_center_lat)
564      max_add = 1
565      do n=1,num_srch_bins
566        if (plat >= bin_lats(1,n) .and. plat <= bin_lats(2,n) .and. &
567            plon >= bin_lons(1,n) .and. plon <= bin_lons(2,n)) then
568          min_add = min(min_add, src_bin_add(1,n))
569          max_add = max(max_add, src_bin_add(2,n))
570        endif
571      end do
572 
573!-----------------------------------------------------------------------
574!
575!     now perform a more detailed search
576!
577!-----------------------------------------------------------------------
578
579      nx = src_grid_dims(1)
580      ny = src_grid_dims(2)
581
582      srch_loop: do srch_add = min_add,max_add
583
584        if (plat <= src_bound_box(2,srch_add) .and.  &
585            plat >= src_bound_box(1,srch_add) .and. &
586            plon <= src_bound_box(4,srch_add) .and.  &
587            plon >= src_bound_box(3,srch_add)) then
588
589          !***
590          !*** we are within bounding box so get really serious
591          !***
592
593          !*** find N,S and NE points to this grid point
594
595          j = (srch_add - 1)/nx +1
596          i = srch_add - (j-1)*nx
597
598          if (i < nx) then
599            ip1 = i + 1
600          else
601            ip1 = 1
602          endif
603
604          if (j < ny) then
605            jp1 = j+1
606          else
607            jp1 = 1
608          endif
609
610          n_add = (jp1 - 1)*nx + i
611          e_add = (j - 1)*nx + ip1
612          ne_add = (jp1 - 1)*nx + ip1
613
614          !***
615          !*** find N,S and NE lat/lon coords and check bounding box
616          !***
617
618          src_lats(1) = src_center_lat(srch_add)
619          src_lats(2) = src_center_lat(e_add)
620          src_lats(3) = src_center_lat(ne_add)
621          src_lats(4) = src_center_lat(n_add)
622
623          src_lons(1) = src_center_lon(srch_add)
624          src_lons(2) = src_center_lon(e_add)
625          src_lons(3) = src_center_lon(ne_add)
626          src_lons(4) = src_center_lon(n_add)
627
628          !***
629          !*** for consistency, we must make sure all lons are in
630          !*** same 2pi interval
631          !***
632
633          vec1_lon = src_lons(1) - plon
634          if (vec1_lon > pi) then
635            src_lons(1) = src_lons(1) - pi2
636          else if (vec1_lon < -pi) then
637            src_lons(1) = src_lons(1) + pi2
638          endif
639          do n=2,4
640            vec1_lon = src_lons(n) - src_lons(1)
641            if (vec1_lon > pi) then
642              src_lons(n) = src_lons(n) - pi2
643            else if (vec1_lon < -pi) then
644              src_lons(n) = src_lons(n) + pi2
645            endif
646          end do
647
648          corner_loop: do n=1,4
649            next_n = MOD(n,4) + 1
650
651            !***
652            !*** here we take the cross product of the vector making
653            !*** up each box side with the vector formed by the vertex
654            !*** and search point.  if all the cross products are
655            !*** same sign, the point is contained in the box.
656            !***
657
658            vec1_lat = src_lats(next_n) - src_lats(n)
659            vec1_lon = src_lons(next_n) - src_lons(n)
660            vec2_lat = plat - src_lats(n)
661            vec2_lon = plon - src_lons(n)
662
663            !***
664            !*** check for 0,2pi crossings
665            !***
666
667            if (vec1_lon >  three*pih) then
668              vec1_lon = vec1_lon - pi2
669            else if (vec1_lon < -three*pih) then
670              vec1_lon = vec1_lon + pi2
671            endif
672            if (vec2_lon >  three*pih) then
673              vec2_lon = vec2_lon - pi2
674            else if (vec2_lon < -three*pih) then
675              vec2_lon = vec2_lon + pi2
676            endif
677
678            cross_product = vec1_lon*vec2_lat - vec2_lon*vec1_lat
679
680            !***
681            !*** if cross product is less than zero, this cell
682            !*** doesn't work
683            !***
684
685            if (n==1) cross_product_last = cross_product
686            if (cross_product*cross_product_last < zero) then
687              exit corner_loop
688            else
689              cross_product_last = cross_product
690            endif
691
692          end do corner_loop
693
694          !***
695          !*** if cross products all positive, we found the location
696          !***
697
698          if (n > 4) then
699            src_add(1) = srch_add
700            src_add(2) = e_add
701            src_add(3) = ne_add
702            src_add(4) = n_add
703
704            return
705          endif
706
707          !***
708          !*** otherwise move on to next cell
709          !***
710
711        endif !bounding box check
712      end do srch_loop
713
714      !***
715      !*** if no cell found, point is likely either in a box that
716      !*** straddles either pole or is outside the grid.  fall back
717      !*** to a distance-weighted average of the four closest
718      !*** points.  go ahead and compute weights here, but store
719      !*** in src_lats and return -add to prevent the parent
720      !*** routine from computing bilinear weights
721      !***
722
723      coslat_dst = cos(plat)
724      sinlat_dst = sin(plat)
725      coslon_dst = cos(plon)
726      sinlon_dst = sin(plon)
727
728      dist_min = bignum
729      src_lats = bignum
730      do srch_add = min_add,max_add
731        distance = acos(coslat_dst*cos(src_center_lat(srch_add))* &
732                       (coslon_dst*cos(src_center_lon(srch_add)) + &
733                        sinlon_dst*sin(src_center_lon(srch_add)))+ &
734                        sinlat_dst*sin(src_center_lat(srch_add)))
735
736        if (distance < dist_min) then
737          sort_loop: do n=1,4
738            if (distance < src_lats(n)) then
739              do i=4,n+1,-1
740                src_add (i) = src_add (i-1)
741                src_lats(i) = src_lats(i-1)
742              end do
743              src_add (n) = -srch_add
744              src_lats(n) = distance
745              dist_min = src_lats(4)
746              exit sort_loop
747            endif
748          end do sort_loop
749        endif
750      end do
751
752      src_lons = one/(src_lats + tiny)
753      distance = sum(src_lons)
754      src_lats = src_lons/distance
755
756!-----------------------------------------------------------------------
757
758      end subroutine grid_search_bicub 
759
760!***********************************************************************
761
762      subroutine store_link_bicub(dst_add, src_add, weights, nmap)
763
764!-----------------------------------------------------------------------
765!
766!     this routine stores the address and weight for four links
767!     associated with one destination point in the appropriate address
768!     and weight arrays and resizes those arrays if necessary.
769!
770!-----------------------------------------------------------------------
771
772!-----------------------------------------------------------------------
773!
774!     input variables
775!
776!-----------------------------------------------------------------------
777
778      integer (kind=int_kind), intent(in) :: &
779              dst_add,   & ! address on destination grid
780              nmap      ! identifies which direction for mapping
781
782      integer (kind=int_kind), dimension(4), intent(in) :: &
783              src_add   ! addresses on source grid
784
785      real (kind=dbl_kind), dimension(4,4), intent(in) :: &
786              weights ! array of remapping weights for these links
787
788!-----------------------------------------------------------------------
789!
790!     local variables
791!
792!-----------------------------------------------------------------------
793
794      integer (kind=int_kind) :: n,  & ! dummy index
795             num_links_old          ! placeholder for old link number
796
797!-----------------------------------------------------------------------
798!
799!     increment number of links and check to see if remap arrays need
800!     to be increased to accomodate the new link.  then store the
801!     link.
802!
803!-----------------------------------------------------------------------
804
805      select case (nmap)
806      case(1)
807
808        num_links_old  = num_links_map1
809        num_links_map1 = num_links_old + 4
810
811        if (num_links_map1 > max_links_map1)  &
812           call resize_remap_vars(1,resize_increment)
813
814        do n=1,4
815          grid1_add_map1(num_links_old+n) = src_add(n)
816          grid2_add_map1(num_links_old+n) = dst_add
817          wts_map1    (:,num_links_old+n) = weights(:,n)
818        end do
819
820      case(2)
821
822        num_links_old  = num_links_map2
823        num_links_map2 = num_links_old + 4
824
825        if (num_links_map2 > max_links_map2)  &
826           call resize_remap_vars(2,resize_increment)
827
828        do n=1,4
829          grid1_add_map2(num_links_old+n) = dst_add
830          grid2_add_map2(num_links_old+n) = src_add(n)
831          wts_map2    (:,num_links_old+n) = weights(:,n)
832        end do
833
834      end select
835
836!-----------------------------------------------------------------------
837
838      end subroutine store_link_bicub
839
840!***********************************************************************
841
842      end module remap_bicubic
843
844!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Note: See TracBrowser for help on using the repository browser.