source: branches/UKMO/r6232_tracer_advection/NEMOGCM/TOOLS/WEIGHTS/src/remap_distwgt.f90 @ 9295

Last change on this file since 9295 was 9295, checked in by jcastill, 3 years ago

Remove svn keywords

File size: 15.6 KB
Line 
1!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2!
3!     this module contains necessary routines for performing an
4!     interpolation using a distance-weighted average of n nearest
5!     neighbors.
6!
7!-----------------------------------------------------------------------
8!
9!     CVS:$Id$
10!
11!     Copyright (c) 1997, 1998 the Regents of the University of
12!       California.
13!
14!     This software and ancillary information (herein called software)
15!     called SCRIP is made available under the terms described here. 
16!     The software has been approved for release with associated
17!     LA-CC Number 98-45.
18!
19!     Unless otherwise indicated, this software has been authored
20!     by an employee or employees of the University of California,
21!     operator of the Los Alamos National Laboratory under Contract
22!     No. W-7405-ENG-36 with the U.S. Department of Energy.  The U.S.
23!     Government has rights to use, reproduce, and distribute this
24!     software.  The public may copy and use this software without
25!     charge, provided that this Notice and any statement of authorship
26!     are reproduced on all copies.  Neither the Government nor the
27!     University makes any warranty, express or implied, or assumes
28!     any liability or responsibility for the use of this software.
29!
30!     If software is modified to produce derivative works, such modified
31!     software should be clearly marked, so as not to confuse it with
32!     the version available from Los Alamos National Laboratory.
33!
34!***********************************************************************
35
36      module remap_distance_weight
37
38!-----------------------------------------------------------------------
39
40      use kinds_mod     ! defines common data types
41      use constants     ! defines common constants
42      use grids         ! module containing grid info
43      use remap_vars    ! module containing remap info
44
45      implicit none
46
47!-----------------------------------------------------------------------
48!
49!     module variables
50!
51!-----------------------------------------------------------------------
52
53      integer (kind=int_kind), parameter ::  &
54           num_neighbors=4  ! num nearest neighbors to interpolate from
55
56      real (kind=dbl_kind), dimension(:), allocatable, save :: &
57           coslat, sinlat,  & ! cosine, sine of grid lats (for distance)
58           coslon, sinlon,  & ! cosine, sine of grid lons (for distance)
59           wgtstmp         ! an array to hold the link weight
60
61!***********************************************************************
62
63      contains
64
65!***********************************************************************
66
67      subroutine remap_distwgt
68
69!-----------------------------------------------------------------------
70!
71!     this routine computes the inverse-distance weights for a
72!     nearest-neighbor interpolation.
73!
74!-----------------------------------------------------------------------
75
76!-----------------------------------------------------------------------
77!
78!     local variables
79!
80!-----------------------------------------------------------------------
81
82      logical (kind=log_kind), dimension(num_neighbors) :: &
83           nbr_mask        ! mask at nearest neighbors
84
85      integer (kind=int_kind) :: n, &
86           dst_add,         & ! destination address
87           nmap            ! index of current map being computed
88
89      integer (kind=int_kind), dimension(num_neighbors) :: &
90           nbr_add         ! source address at nearest neighbors
91
92      real (kind=dbl_kind), dimension(num_neighbors) :: &
93           nbr_dist        ! angular distance four nearest neighbors
94
95      real (kind=dbl_kind) :: &
96           coslat_dst,      & ! cos(lat) of destination grid point
97           coslon_dst,      & ! cos(lon) of destination grid point
98           sinlat_dst,      & ! sin(lat) of destination grid point
99           sinlon_dst,      & ! sin(lon) of destination grid point
100           dist_tot        ! sum of neighbor distances (for normalizing)
101
102!-----------------------------------------------------------------------
103!
104!     compute mappings from grid1 to grid2
105!
106!-----------------------------------------------------------------------
107
108      nmap = 1
109
110      !***
111      !*** allocate wgtstmp to be consistent with store_link interface
112      !***
113
114      allocate (wgtstmp(num_wts))
115
116      !***
117      !*** compute cos, sin of lat/lon on source grid for distance
118      !*** calculations
119      !***
120
121      allocate (coslat(grid1_size), coslon(grid1_size), &
122                sinlat(grid1_size), sinlon(grid1_size))
123
124      coslat = cos(grid1_center_lat)
125      coslon = cos(grid1_center_lon)
126      sinlat = sin(grid1_center_lat)
127      sinlon = sin(grid1_center_lon)
128
129      !***
130      !*** loop over destination grid
131      !***
132
133      grid_loop1: do dst_add = 1, grid2_size
134
135        if (.not. grid2_mask(dst_add)) cycle grid_loop1
136
137        coslat_dst = cos(grid2_center_lat(dst_add))
138        coslon_dst = cos(grid2_center_lon(dst_add))
139        sinlat_dst = sin(grid2_center_lat(dst_add))
140        sinlon_dst = sin(grid2_center_lon(dst_add))
141
142        !***
143        !*** find nearest grid points on source grid and
144        !*** distances to each point
145        !***
146
147        call grid_search_nbr(nbr_add, nbr_dist,  &
148                             grid2_center_lat(dst_add), &
149                             grid2_center_lon(dst_add), &
150                             coslat_dst, coslon_dst,  &
151                             sinlat_dst, sinlon_dst, &
152                             bin_addr1, bin_addr2)
153
154        !***
155        !*** compute weights based on inverse distance
156        !*** if mask is false, eliminate those points
157        !***
158
159        dist_tot = zero
160        do n=1,num_neighbors
161          if (grid1_mask(nbr_add(n))) then
162            nbr_dist(n) = one/nbr_dist(n)
163            dist_tot = dist_tot + nbr_dist(n)
164            nbr_mask(n) = .true.
165          else
166            nbr_mask(n) = .false.
167          endif
168        end do
169
170        !***
171        !*** normalize weights and store the link
172        !***
173
174        do n=1,num_neighbors
175          if (nbr_mask(n)) then
176            wgtstmp(1) = nbr_dist(n)/dist_tot
177            call store_link_nbr(nbr_add(n), dst_add, wgtstmp, nmap)
178            grid2_frac(dst_add) = one
179          endif
180        end do
181
182      end do grid_loop1
183
184      deallocate (coslat, coslon, sinlat, sinlon)
185
186!-----------------------------------------------------------------------
187!
188!     compute mappings from grid2 to grid1 if necessary
189!
190!-----------------------------------------------------------------------
191
192      if (num_maps > 1) then
193
194      nmap = 2
195
196      !***
197      !*** compute cos, sin of lat/lon on source grid for distance
198      !*** calculations
199      !***
200
201      allocate (coslat(grid2_size), coslon(grid2_size), &
202                sinlat(grid2_size), sinlon(grid2_size))
203
204      coslat = cos(grid2_center_lat)
205      coslon = cos(grid2_center_lon)
206      sinlat = sin(grid2_center_lat)
207      sinlon = sin(grid2_center_lon)
208
209      !***
210      !*** loop over destination grid
211      !***
212
213      grid_loop2: do dst_add = 1, grid1_size
214
215        if (.not. grid1_mask(dst_add)) cycle grid_loop2
216
217        coslat_dst = cos(grid1_center_lat(dst_add))
218        coslon_dst = cos(grid1_center_lon(dst_add))
219        sinlat_dst = sin(grid1_center_lat(dst_add))
220        sinlon_dst = sin(grid1_center_lon(dst_add))
221
222        !***
223        !*** find four nearest grid points on source grid and
224        !*** distances to each point
225        !***
226
227        call grid_search_nbr(nbr_add, nbr_dist, &
228                             grid1_center_lat(dst_add), &
229                             grid1_center_lon(dst_add), &
230                             coslat_dst, coslon_dst,  &
231                             sinlat_dst, sinlon_dst, &
232                             bin_addr2, bin_addr1)
233
234        !***
235        !*** compute weights based on inverse distance
236        !*** if mask is false, eliminate those points
237        !***
238
239        dist_tot = zero
240        do n=1,num_neighbors
241          if (grid2_mask(nbr_add(n))) then
242            nbr_dist(n) = one/nbr_dist(n)
243            dist_tot = dist_tot + nbr_dist(n)
244            nbr_mask(n) = .true.
245          else
246            nbr_mask(n) = .false.
247          endif
248        end do
249
250        !***
251        !*** normalize weights and store the link
252        !***
253
254        do n=1,num_neighbors
255          if (nbr_mask(n)) then
256            wgtstmp(1) = nbr_dist(n)/dist_tot
257            call store_link_nbr(dst_add, nbr_add(n), wgtstmp, nmap)
258            grid1_frac(dst_add) = one
259          endif
260        end do
261
262      end do grid_loop2
263
264      deallocate (coslat, coslon, sinlat, sinlon)
265
266      endif
267
268      deallocate(wgtstmp)
269
270!-----------------------------------------------------------------------
271
272      end subroutine remap_distwgt
273
274!***********************************************************************
275
276      subroutine grid_search_nbr(nbr_add, nbr_dist, plat, plon,  &
277                     coslat_dst, coslon_dst, sinlat_dst, sinlon_dst, &
278                     src_bin_add, dst_bin_add)
279
280!-----------------------------------------------------------------------
281!
282!     this routine finds the closest num_neighbor points to a search
283!     point and computes a distance to each of the neighbors.
284!
285!-----------------------------------------------------------------------
286
287!-----------------------------------------------------------------------
288!
289!     output variables
290!
291!-----------------------------------------------------------------------
292
293      integer (kind=int_kind), dimension(num_neighbors), intent(out) :: &
294              nbr_add  ! address of each of the closest points
295
296      real (kind=dbl_kind), dimension(num_neighbors), intent(out) :: &
297              nbr_dist ! distance to each of the closest points
298
299!-----------------------------------------------------------------------
300!
301!     input variables
302!
303!-----------------------------------------------------------------------
304
305      integer (kind=int_kind), dimension(:,:), intent(in) :: &
306              src_bin_add,  & ! search bins for restricting search
307              dst_bin_add   
308
309      real (kind=dbl_kind), intent(in) :: &
310              plat,          & ! latitude  of the search point
311              plon,          & ! longitude of the search point
312              coslat_dst,    & ! cos(lat)  of the search point
313              coslon_dst,    & ! cos(lon)  of the search point
314              sinlat_dst,    & ! sin(lat)  of the search point
315              sinlon_dst    ! sin(lon)  of the search point
316
317!-----------------------------------------------------------------------
318!
319!     local variables
320!
321!-----------------------------------------------------------------------
322
323      integer (kind=int_kind) :: n, nmax, nadd, nchk,  & ! dummy indices
324              min_add, max_add, nm1, np1, i, j, ip1, im1, jp1, jm1
325
326      real (kind=dbl_kind) :: &
327              distance      ! angular distance
328
329!-----------------------------------------------------------------------
330!
331!     loop over source grid and find nearest neighbors
332!
333!-----------------------------------------------------------------------
334
335      !***
336      !*** restrict the search using search bins
337      !*** expand the bins to catch neighbors
338      !***
339
340      select case (restrict_type)
341      case('latitude')
342
343        do n=1,num_srch_bins
344          if (plat >= bin_lats(1,n) .and. plat <= bin_lats(2,n)) then
345            min_add = src_bin_add(1,n)
346            max_add = src_bin_add(2,n)
347
348            nm1 = max(n-1,1)
349            np1 = min(n+1,num_srch_bins)
350
351            min_add = min(min_add,src_bin_add(1,nm1))
352            max_add = max(max_add,src_bin_add(2,nm1))
353            min_add = min(min_add,src_bin_add(1,np1))
354            max_add = max(max_add,src_bin_add(2,np1))
355          endif
356        end do
357
358      case('latlon')
359
360        n = 0
361        nmax = nint(sqrt(real(num_srch_bins)))
362        do j=1,nmax
363        jp1 = min(j+1,nmax)
364        jm1 = max(j-1,1)
365        do i=1,nmax
366          ip1 = min(i+1,nmax)
367          im1 = max(i-1,1)
368
369          n = n+1
370          if (plat >= bin_lats(1,n) .and. plat <= bin_lats(2,n) .and. &
371              plon >= bin_lons(1,n) .and. plon <= bin_lons(3,n)) then
372            min_add = src_bin_add(1,n)
373            max_add = src_bin_add(2,n)
374
375            nm1 = (jm1-1)*nmax + im1
376            np1 = (jp1-1)*nmax + ip1
377            nm1 = max(nm1,1)
378            np1 = min(np1,num_srch_bins)
379
380            min_add = min(min_add,src_bin_add(1,nm1))
381            max_add = max(max_add,src_bin_add(2,nm1))
382            min_add = min(min_add,src_bin_add(1,np1))
383            max_add = max(max_add,src_bin_add(2,np1))
384          endif
385        end do
386        end do
387
388      end select
389
390      !***
391      !*** initialize distance and address arrays
392      !***
393
394      nbr_add = 0
395      nbr_dist = bignum
396
397      do nadd=min_add,max_add
398
399        !***
400        !*** find distance to this point
401        !***
402
403        distance = acos(sinlat_dst*sinlat(nadd) + &
404                        coslat_dst*coslat(nadd)* &
405                       (coslon_dst*coslon(nadd) + &
406                        sinlon_dst*sinlon(nadd)) )
407
408        !***
409        !*** store the address and distance if this is one of the
410        !*** smallest four so far
411        !***
412
413        check_loop: do nchk=1,num_neighbors
414          if (distance .lt. nbr_dist(nchk)) then
415            do n=num_neighbors,nchk+1,-1
416              nbr_add(n) = nbr_add(n-1)
417              nbr_dist(n) = nbr_dist(n-1)
418            end do
419            nbr_add(nchk) = nadd
420            nbr_dist(nchk) = distance
421            exit check_loop
422          endif
423        end do check_loop
424
425      end do
426
427!-----------------------------------------------------------------------
428
429      end subroutine grid_search_nbr 
430
431!***********************************************************************
432
433      subroutine store_link_nbr(add1, add2, weights, nmap)
434
435!-----------------------------------------------------------------------
436!
437!     this routine stores the address and weight for this link in
438!     the appropriate address and weight arrays and resizes those
439!     arrays if necessary.
440!
441!-----------------------------------------------------------------------
442
443!-----------------------------------------------------------------------
444!
445!     input variables
446!
447!-----------------------------------------------------------------------
448
449      integer (kind=int_kind), intent(in) :: &
450              add1,   & ! address on grid1
451              add2,   & ! address on grid2
452              nmap   ! identifies which direction for mapping
453
454      real (kind=dbl_kind), dimension(:), intent(in) :: &
455              weights ! array of remapping weights for this link
456
457!-----------------------------------------------------------------------
458!
459!     increment number of links and check to see if remap arrays need
460!     to be increased to accomodate the new link.  then store the
461!     link.
462!
463!-----------------------------------------------------------------------
464
465      select case (nmap)
466      case(1)
467
468        num_links_map1  = num_links_map1 + 1
469
470        if (num_links_map1 > max_links_map1)  &
471           call resize_remap_vars(1,resize_increment)
472
473        grid1_add_map1(num_links_map1) = add1
474        grid2_add_map1(num_links_map1) = add2
475        wts_map1    (:,num_links_map1) = weights
476
477      case(2)
478
479        num_links_map2  = num_links_map2 + 1
480
481        if (num_links_map2 > max_links_map2)  &
482           call resize_remap_vars(2,resize_increment)
483
484        grid1_add_map2(num_links_map2) = add1
485        grid2_add_map2(num_links_map2) = add2
486        wts_map2    (:,num_links_map2) = weights
487
488      end select
489
490!-----------------------------------------------------------------------
491
492      end subroutine store_link_nbr
493
494!***********************************************************************
495
496      end module remap_distance_weight
497
498!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Note: See TracBrowser for help on using the repository browser.