source: branches/2015/dev_r5021_UKMO1_CICE_coupling/NEMOGCM/TOOLS/WEIGHTS/src/remap_bilinear.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: 24.5 KB
Line 
1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2!
3!     this module contains necessary routines for performing an
4!     bilinear 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_bilinear
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_bilin
61
62!-----------------------------------------------------------------------
63!
64!     this routine computes the weights for a bilinear 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           wgts            ! bilinear weights for four corners
86
87      real (kind=dbl_kind) :: &
88           plat, plon,        & ! lat/lon coords of destination point
89           iguess, jguess,    & ! current guess for bilinear coordinate
90           thguess, phguess,  & ! current guess for lat/lon coordinate
91           deli, delj,        & ! corrections to i,j
92           dth1, dth2, dth3,  & ! some latitude  differences
93           dph1, dph2, dph3,  & ! some longitude differences
94           dthp, dphp,        & ! difference between point and sw corner
95           mat1, mat2, mat3, mat4,  & ! matrix elements
96           determinant,       & ! matrix determinant
97           sum_wgts          ! sum of weights for normalization
98
99!-----------------------------------------------------------------------
100!
101!     compute mappings from grid1 to grid2
102!
103!-----------------------------------------------------------------------
104
105      nmap = 1
106      if (grid1_rank /= 2) then
107        stop 'Can not do bilinear interpolation when grid_rank /= 2'
108      endif
109
110      !***
111      !*** loop over destination grid
112      !***
113
114      grid_loop1: do dst_add = 1, grid2_size
115
116        if (.not. grid2_mask(dst_add)) cycle grid_loop1
117
118        plat = grid2_center_lat(dst_add)
119        plon = grid2_center_lon(dst_add)
120
121        !***
122        !*** find nearest square of grid points on source grid
123        !***
124
125        call grid_search_bilin(src_add, src_lats, src_lons,  &
126                               plat, plon, grid1_dims, &
127                               grid1_center_lat, grid1_center_lon, &
128                               grid1_bound_box, bin_addr1, bin_addr2)
129
130        !***
131        !*** check to see if points are land points
132        !***
133
134        if (src_add(1) > 0) then
135          do n=1,4
136            if (.not. grid1_mask(src_add(n))) src_add(1) = 0
137          end do
138        endif
139
140        !***
141        !*** if point found, find local i,j coordinates for weights
142        !***
143
144        if (src_add(1) > 0) then
145
146          grid2_frac(dst_add) = one
147
148          !***
149          !*** iterate to find i,j for bilinear approximation
150          !***
151
152          dth1 = src_lats(2) - src_lats(1)
153          dth2 = src_lats(4) - src_lats(1)
154          dth3 = src_lats(3) - src_lats(2) - dth2
155
156          dph1 = src_lons(2) - src_lons(1)
157          dph2 = src_lons(4) - src_lons(1)
158          dph3 = src_lons(3) - src_lons(2)
159
160          if (dph1 >  three*pih) dph1 = dph1 - pi2
161          if (dph2 >  three*pih) dph2 = dph2 - pi2
162          if (dph3 >  three*pih) dph3 = dph3 - pi2
163          if (dph1 < -three*pih) dph1 = dph1 + pi2
164          if (dph2 < -three*pih) dph2 = dph2 + pi2
165          if (dph3 < -three*pih) dph3 = dph3 + pi2
166
167          dph3 = dph3 - dph2
168
169          iguess = half
170          jguess = half
171
172          iter_loop1: do iter=1,max_iter
173
174            dthp = plat - src_lats(1) - dth1*iguess - &
175                          dth2*jguess - dth3*iguess*jguess
176            dphp = plon - src_lons(1)
177
178            if (dphp >  three*pih) dphp = dphp - pi2
179            if (dphp < -three*pih) dphp = dphp + pi2
180
181            dphp = dphp - dph1*iguess - dph2*jguess -  &
182                          dph3*iguess*jguess
183
184            mat1 = dth1 + dth3*jguess
185            mat2 = dth2 + dth3*iguess
186            mat3 = dph1 + dph3*jguess
187            mat4 = dph2 + dph3*iguess
188
189            determinant = mat1*mat4 - mat2*mat3
190
191            deli = (dthp*mat4 - mat2*dphp)/determinant
192            delj = (mat1*dphp - dthp*mat3)/determinant
193
194            if (abs(deli) < converge .and.  &
195                abs(delj) < converge) exit iter_loop1
196
197            iguess = iguess + deli
198            jguess = jguess + delj
199
200          end do iter_loop1
201
202          if (iter <= max_iter) then
203
204            !***
205            !*** successfully found i,j - compute weights
206            !***
207
208            wgts(1) = (one-iguess)*(one-jguess)
209            wgts(2) = iguess*(one-jguess)
210            wgts(3) = iguess*jguess
211            wgts(4) = (one-iguess)*jguess
212
213            call store_link_bilin(dst_add, src_add, wgts, nmap)
214
215          else
216            print *,'Point coords: ',plat,plon
217            print *,'Dest grid lats: ',src_lats
218            print *,'Dest grid lons: ',src_lons
219            print *,'Dest grid addresses: ',src_add
220            print *,'Current i,j : ',iguess, jguess
221            stop 'Iteration for i,j exceed max iteration count'
222          endif
223
224        !***
225        !*** search for bilinear failed - use a distance-weighted
226        !*** average instead (this is typically near the pole)
227        !***
228
229        else if (src_add(1) < 0) then
230
231          src_add = abs(src_add)
232          icount = 0
233          do n=1,4
234            if (grid1_mask(src_add(n))) then
235              icount = icount + 1
236            else
237              src_lats(n) = zero
238            endif
239          end do
240
241          if (icount > 0) then
242            !*** renormalize weights
243
244            sum_wgts = sum(src_lats)
245            wgts(1) = src_lats(1)/sum_wgts
246            wgts(2) = src_lats(2)/sum_wgts
247            wgts(3) = src_lats(3)/sum_wgts
248            wgts(4) = src_lats(4)/sum_wgts
249
250            grid2_frac(dst_add) = one
251            call store_link_bilin(dst_add, src_add, wgts, nmap)
252          endif
253
254        endif
255      end do grid_loop1
256
257!-----------------------------------------------------------------------
258!
259!     compute mappings from grid2 to grid1 if necessary
260!
261!-----------------------------------------------------------------------
262
263      if (num_maps > 1) then
264
265      nmap = 2
266      if (grid2_rank /= 2) then
267        stop 'Can not do bilinear interpolation when grid_rank /= 2'
268      endif
269
270      !***
271      !*** loop over destination grid
272      !***
273
274      grid_loop2: do dst_add = 1, grid1_size
275
276        if (.not. grid1_mask(dst_add)) cycle grid_loop2
277
278        plat = grid1_center_lat(dst_add)
279        plon = grid1_center_lon(dst_add)
280
281        !***
282        !*** find nearest square of grid points on source grid
283        !***
284
285        call grid_search_bilin(src_add, src_lats, src_lons,  &
286                               plat, plon, grid2_dims, &
287                               grid2_center_lat, grid2_center_lon, &
288                               grid2_bound_box, bin_addr2, bin_addr1)
289
290        !***
291        !*** check to see if points are land points
292        !***
293
294        if (src_add(1) > 0) then
295          do n=1,4
296            if (.not. grid2_mask(src_add(n))) src_add(1) = 0
297          end do
298        endif
299
300        !***
301        !*** if point found, find i,j coordinates for weights
302        !***
303
304        if (src_add(1) > 0) then
305
306          grid1_frac(dst_add) = one
307
308          !***
309          !*** iterate to find i,j for bilinear approximation
310          !***
311
312          dth1 = src_lats(2) - src_lats(1)
313          dth2 = src_lats(4) - src_lats(1)
314          dth3 = src_lats(3) - src_lats(2) - dth2
315
316          dph1 = src_lons(2) - src_lons(1)
317          dph2 = src_lons(4) - src_lons(1)
318          dph3 = src_lons(3) - src_lons(2)
319
320          if (dph1 >  pi) dph1 = dph1 - pi2
321          if (dph2 >  pi) dph2 = dph2 - pi2
322          if (dph3 >  pi) dph3 = dph3 - pi2
323          if (dph1 < -pi) dph1 = dph1 + pi2
324          if (dph2 < -pi) dph2 = dph2 + pi2
325          if (dph3 < -pi) dph3 = dph3 + pi2
326
327          dph3 = dph3 - dph2
328
329          iguess = zero
330          jguess = zero
331
332          iter_loop2: do iter=1,max_iter
333
334            dthp = plat - src_lats(1) - dth1*iguess - &
335                          dth2*jguess - dth3*iguess*jguess
336            dphp = plon - src_lons(1)
337
338            if (dphp >  pi) dphp = dphp - pi2
339            if (dphp < -pi) dphp = dphp + pi2
340
341            dphp = dphp - dph1*iguess - dph2*jguess -  &
342                          dph3*iguess*jguess
343
344            mat1 = dth1 + dth3*jguess
345            mat2 = dth2 + dth3*iguess
346            mat3 = dph1 + dph3*jguess
347            mat4 = dph2 + dph3*iguess
348
349            determinant = mat1*mat4 - mat2*mat3
350
351            deli = (dthp*mat4 - mat2*dphp)/determinant
352            delj = (mat1*dphp - dthp*mat3)/determinant
353
354            if (abs(deli) < converge .and.  &
355                abs(delj) < converge) exit iter_loop2
356
357            iguess = iguess + deli
358            jguess = jguess + delj
359
360          end do iter_loop2
361
362          if (iter <= max_iter) then
363
364            !***
365            !*** successfully found i,j - compute weights
366            !***
367
368            wgts(1) = (one-iguess)*(one-jguess)
369            wgts(2) = iguess*(one-jguess)
370            wgts(3) = iguess*jguess
371            wgts(4) = (one-iguess)*jguess
372
373            call store_link_bilin(dst_add, src_add, wgts, nmap)
374
375          else
376            print *,'Point coords: ',plat,plon
377            print *,'Dest grid lats: ',src_lats
378            print *,'Dest grid lons: ',src_lons
379            print *,'Dest grid addresses: ',src_add
380            print *,'Current i,j : ',iguess, jguess
381            stop 'Iteration for i,j exceed max iteration count'
382          endif
383
384        !***
385        !*** search for bilinear failed - us a distance-weighted
386        !*** average instead
387        !***
388
389        else if (src_add(1) < 0) then
390
391          src_add = abs(src_add)
392          icount = 0
393          do n=1,4
394            if (grid2_mask(src_add(n))) then
395              icount = icount + 1
396            else
397              src_lats(n) = zero
398            endif
399          end do
400
401          if (icount > 0) then
402            !*** renormalize weights
403
404            sum_wgts = sum(src_lats)
405            wgts(1) = src_lats(1)/sum_wgts
406            wgts(2) = src_lats(2)/sum_wgts
407            wgts(3) = src_lats(3)/sum_wgts
408            wgts(4) = src_lats(4)/sum_wgts
409
410            grid1_frac(dst_add) = one
411            call store_link_bilin(dst_add, src_add, wgts, nmap)
412          endif
413
414        endif
415      end do grid_loop2
416
417      endif ! nmap=2
418
419!-----------------------------------------------------------------------
420
421      end subroutine remap_bilin
422
423!***********************************************************************
424
425      subroutine grid_search_bilin(src_add, src_lats, src_lons,  &
426                                   plat, plon, src_grid_dims, &
427                                   src_center_lat, src_center_lon, &
428                                   src_grid_bound_box, &
429                                   src_bin_add, dst_bin_add)
430
431!-----------------------------------------------------------------------
432!
433!     this routine finds the location of the search point plat, plon
434!     in the source grid and returns the corners needed for a bilinear
435!     interpolation.
436!
437!-----------------------------------------------------------------------
438
439!-----------------------------------------------------------------------
440!
441!     output variables
442!
443!-----------------------------------------------------------------------
444
445      integer (kind=int_kind), dimension(4), intent(out) :: &
446              src_add  ! address of each corner point enclosing P
447
448      real (kind=dbl_kind), dimension(4), intent(out) :: &
449              src_lats,  & ! latitudes  of the four corner points
450              src_lons  ! longitudes of the four corner points
451
452!-----------------------------------------------------------------------
453!
454!     input variables
455!
456!-----------------------------------------------------------------------
457
458      real (kind=dbl_kind), intent(in) :: &
459              plat,    & ! latitude  of the search point
460              plon    ! longitude of the search point
461
462      integer (kind=int_kind), dimension(2), intent(in) :: &
463              src_grid_dims  ! size of each src grid dimension
464
465      real (kind=dbl_kind), dimension(:), intent(in) :: &
466              src_center_lat,  & ! latitude  of each src grid center
467              src_center_lon  ! longitude of each src grid center
468
469      real (kind=dbl_kind), dimension(:,:), intent(in) :: &
470              src_grid_bound_box ! bound box for source grid
471
472      integer (kind=int_kind), dimension(:,:), intent(in) :: &
473              src_bin_add,     & ! latitude bins for restricting
474              dst_bin_add     ! searches
475
476!-----------------------------------------------------------------------
477!
478!     local variables
479!
480!-----------------------------------------------------------------------
481
482      integer (kind=int_kind) :: n, next_n, srch_add,    & ! dummy indices
483          nx, ny,             & ! dimensions of src grid
484          min_add, max_add,   & ! addresses for restricting search
485          i, j, jp1, ip1, n_add, e_add, ne_add  ! addresses
486
487      real (kind=dbl_kind) ::   & ! vectors for cross-product check
488            vec1_lat, vec1_lon, &
489            vec2_lat, vec2_lon, cross_product, cross_product_last, &
490            coslat_dst, sinlat_dst, coslon_dst, sinlon_dst, &
491            dist_min, distance ! for computing dist-weighted avg
492
493!-----------------------------------------------------------------------
494!
495!     restrict search first using bins
496!
497!-----------------------------------------------------------------------
498
499      src_add = 0
500
501      min_add = size(src_center_lat)
502      max_add = 1
503      do n=1,num_srch_bins
504        if (plat >= bin_lats(1,n) .and. plat <= bin_lats(2,n) .and. &
505            plon >= bin_lons(1,n) .and. plon <= bin_lons(2,n)) then
506          min_add = min(min_add, src_bin_add(1,n))
507          max_add = max(max_add, src_bin_add(2,n))
508        endif
509      end do
510 
511!-----------------------------------------------------------------------
512!
513!     now perform a more detailed search
514!
515!-----------------------------------------------------------------------
516
517      nx = src_grid_dims(1)
518      ny = src_grid_dims(2)
519
520      srch_loop: do srch_add = min_add,max_add
521
522        !*** first check bounding box
523
524        if (plat <= src_grid_bound_box(2,srch_add) .and.  &
525            plat >= src_grid_bound_box(1,srch_add) .and. &
526            plon <= src_grid_bound_box(4,srch_add) .and.  &
527            plon >= src_grid_bound_box(3,srch_add)) then
528
529          !***
530          !*** we are within bounding box so get really serious
531          !***
532
533          !*** determine neighbor addresses
534
535          j = (srch_add - 1)/nx +1
536          i = srch_add - (j-1)*nx
537
538          if (i < nx) then
539            ip1 = i + 1
540          else
541            ip1 = 1
542          endif
543
544          if (j < ny) then
545            jp1 = j+1
546          else
547            jp1 = 1
548          endif
549
550          n_add = (jp1 - 1)*nx + i
551          e_add = (j - 1)*nx + ip1
552          ne_add = (jp1 - 1)*nx + ip1
553
554          src_lats(1) = src_center_lat(srch_add)
555          src_lats(2) = src_center_lat(e_add)
556          src_lats(3) = src_center_lat(ne_add)
557          src_lats(4) = src_center_lat(n_add)
558
559          src_lons(1) = src_center_lon(srch_add)
560          src_lons(2) = src_center_lon(e_add)
561          src_lons(3) = src_center_lon(ne_add)
562          src_lons(4) = src_center_lon(n_add)
563
564          !***
565          !*** for consistency, we must make sure all lons are in
566          !*** same 2pi interval
567          !***
568
569          vec1_lon = src_lons(1) - plon
570          if (vec1_lon >  pi) then
571            src_lons(1) = src_lons(1) - pi2
572          else if (vec1_lon < -pi) then
573            src_lons(1) = src_lons(1) + pi2
574          endif
575          do n=2,4
576            vec1_lon = src_lons(n) - src_lons(1)
577            if (vec1_lon >  pi) then
578              src_lons(n) = src_lons(n) - pi2
579            else if (vec1_lon < -pi) then
580              src_lons(n) = src_lons(n) + pi2
581            endif
582          end do
583
584          corner_loop: do n=1,4
585            next_n = MOD(n,4) + 1
586
587            !***
588            !*** here we take the cross product of the vector making
589            !*** up each box side with the vector formed by the vertex
590            !*** and search point.  if all the cross products are
591            !*** positive, the point is contained in the box.
592            !***
593
594            vec1_lat = src_lats(next_n) - src_lats(n)
595            vec1_lon = src_lons(next_n) - src_lons(n)
596            vec2_lat = plat - src_lats(n)
597            vec2_lon = plon - src_lons(n)
598
599            !***
600            !*** check for 0,2pi crossings
601            !***
602
603            if (vec1_lon >  three*pih) then
604              vec1_lon = vec1_lon - pi2
605            else if (vec1_lon < -three*pih) then
606              vec1_lon = vec1_lon + pi2
607            endif
608            if (vec2_lon >  three*pih) then
609              vec2_lon = vec2_lon - pi2
610            else if (vec2_lon < -three*pih) then
611              vec2_lon = vec2_lon + pi2
612            endif
613
614            cross_product = vec1_lon*vec2_lat - vec2_lon*vec1_lat
615
616            !***
617            !*** if cross product is less than zero, this cell
618            !*** doesn't work
619            !***
620
621            if (n == 1) cross_product_last = cross_product
622            if (cross_product*cross_product_last < zero)  &
623                exit corner_loop
624            cross_product_last = cross_product
625
626          end do corner_loop
627
628          !***
629          !*** if cross products all same sign, we found the location
630          !***
631
632          if (n > 4) then
633            src_add(1) = srch_add
634            src_add(2) = e_add
635            src_add(3) = ne_add
636            src_add(4) = n_add
637
638            return
639          endif
640
641          !***
642          !*** otherwise move on to next cell
643          !***
644
645        endif !bounding box check
646      end do srch_loop
647
648      !***
649      !*** if no cell found, point is likely either in a box that
650      !*** straddles either pole or is outside the grid.  fall back
651      !*** to a distance-weighted average of the four closest
652      !*** points.  go ahead and compute weights here, but store
653      !*** in src_lats and return -add to prevent the parent
654      !*** routine from computing bilinear weights
655      !***
656
657      !print *,'Could not find location for ',plat,plon
658      !print *,'Using nearest-neighbor average for this point'
659
660      coslat_dst = cos(plat)
661      sinlat_dst = sin(plat)
662      coslon_dst = cos(plon)
663      sinlon_dst = sin(plon)
664
665      dist_min = bignum
666      src_lats = bignum
667      do srch_add = min_add,max_add
668        distance = acos(coslat_dst*cos(src_center_lat(srch_add))* &
669                       (coslon_dst*cos(src_center_lon(srch_add)) + &
670                        sinlon_dst*sin(src_center_lon(srch_add)))+ &
671                        sinlat_dst*sin(src_center_lat(srch_add)))
672
673        if (distance < dist_min) then
674          sort_loop: do n=1,4
675            if (distance < src_lats(n)) then
676              do i=4,n+1,-1
677                src_add (i) = src_add (i-1)
678                src_lats(i) = src_lats(i-1)
679              end do
680              src_add (n) = -srch_add
681              src_lats(n) = distance
682              dist_min = src_lats(4)
683              exit sort_loop
684            endif
685          end do sort_loop
686        endif
687      end do
688
689      src_lons = one/(src_lats + tiny)
690      distance = sum(src_lons)
691      src_lats = src_lons/distance
692
693!-----------------------------------------------------------------------
694
695      end subroutine grid_search_bilin 
696
697!***********************************************************************
698
699      subroutine store_link_bilin(dst_add, src_add, weights, nmap)
700
701!-----------------------------------------------------------------------
702!
703!     this routine stores the address and weight for four links
704!     associated with one destination point in the appropriate address
705!     and weight arrays and resizes those arrays if necessary.
706!
707!-----------------------------------------------------------------------
708
709!-----------------------------------------------------------------------
710!
711!     input variables
712!
713!-----------------------------------------------------------------------
714
715      integer (kind=int_kind), intent(in) :: &
716              dst_add,   & ! address on destination grid
717              nmap      ! identifies which direction for mapping
718
719      integer (kind=int_kind), dimension(4), intent(in) :: &
720              src_add   ! addresses on source grid
721
722      real (kind=dbl_kind), dimension(4), intent(in) :: &
723              weights ! array of remapping weights for these links
724
725!-----------------------------------------------------------------------
726!
727!     local variables
728!
729!-----------------------------------------------------------------------
730
731      integer (kind=int_kind) :: n,  & ! dummy index
732             num_links_old          ! placeholder for old link number
733
734!-----------------------------------------------------------------------
735!
736!     increment number of links and check to see if remap arrays need
737!     to be increased to accomodate the new link.  then store the
738!     link.
739!
740!-----------------------------------------------------------------------
741
742      select case (nmap)
743      case(1)
744
745        num_links_old  = num_links_map1
746        num_links_map1 = num_links_old + 4
747
748        if (num_links_map1 > max_links_map1)  &
749           call resize_remap_vars(1,resize_increment)
750
751        do n=1,4
752          grid1_add_map1(num_links_old+n) = src_add(n)
753          grid2_add_map1(num_links_old+n) = dst_add
754          wts_map1    (1,num_links_old+n) = weights(n)
755        end do
756
757      case(2)
758
759        num_links_old  = num_links_map2
760        num_links_map2 = num_links_old + 4
761
762        if (num_links_map2 > max_links_map2)  &
763           call resize_remap_vars(2,resize_increment)
764
765        do n=1,4
766          grid1_add_map2(num_links_old+n) = dst_add
767          grid2_add_map2(num_links_old+n) = src_add(n)
768          wts_map2    (1,num_links_old+n) = weights(n)
769        end do
770
771      end select
772
773!-----------------------------------------------------------------------
774
775      end subroutine store_link_bilin
776
777!***********************************************************************
778
779      end module remap_bilinear
780
781!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Note: See TracBrowser for help on using the repository browser.