source: CPL/oasis3/trunk/src/lib/scrip/src/remap_distwgt.f @ 1677

Last change on this file since 1677 was 1677, checked in by aclsce, 12 years ago

Imported oasis3 (tag ipslcm5a) from cvs server to svn server (igcmg project).

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