source: CPL/oasis3/trunk/src/lib/scrip/src/remap_gauswgt.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: 21.7 KB
Line 
1C****
2C                    ************************
3C                    *     OASIS MODULE     *
4C                    *     ------------     *
5C                    ************************
6C**** 
7C***********************************************************************
8C     This module belongs to the SCRIP library. It is a modified version
9C     of SCRIP 1.4 remap_ditswgt.f in order to weight a neighbour by
10C     exp[-1/2 * d^^2/sigma^^2] where d is the distance between the source
11C     and the target points and sigma^^2=VAR*dm^^2, where dm is the
12C     average distance between the source grid points and VAR is a value
13C     given by the user.
14C
15C     Additional Modifications:
16C       - restrict types are written in capital letters
17C       - bug line 429: bin_lons(3,n) instead of bin_lons(2,n)
18C
19C     Modified by            D. Declat,  CERFACS              27.06.2002
20C***********************************************************************
21!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
22!
23!     this module contains necessary routines for performing an
24!     interpolation using a distance-weighted average of n nearest
25!     neighbors.
26!
27!-----------------------------------------------------------------------
28!
29!     CVS:$Id: remap_gauswgt.f,v 1.1.1.1 2005/03/23 16:01:12 adm Exp $
30!
31!     Copyright (c) 1997, 1998 the Regents of the University of
32!       California.
33!
34!     This software and ancillary information (herein called software)
35!     called SCRIP is made available under the terms described here. 
36!     The software has been approved for release with associated
37!     LA-CC Number 98-45.
38!
39!     Unless otherwise indicated, this software has been authored
40!     by an employee or employees of the University of California,
41!     operator of the Los Alamos National Laboratory under Contract
42!     No. W-7405-ENG-36 with the U.S. Department of Energy.  The U.S.
43!     Government has rights to use, reproduce, and distribute this
44!     software.  The public may copy and use this software without
45!     charge, provided that this Notice and any statement of authorship
46!     are reproduced on all copies.  Neither the Government nor the
47!     University makes any warranty, express or implied, or assumes
48!     any liability or responsibility for the use of this software.
49!
50!     If software is modified to produce derivative works, such modified
51!     software should be clearly marked, so as not to confuse it with
52!     the version available from Los Alamos National Laboratory.
53!
54!***********************************************************************
55
56      module remap_gaussian_weight
57
58!-----------------------------------------------------------------------
59
60      use kinds_mod     ! defines common data types
61      use constants     ! defines common constants
62      use grids         ! module containing grid info
63      use remap_vars    ! module containing remap info
64
65      implicit none
66
67!-----------------------------------------------------------------------
68!
69!     module variables
70!
71!-----------------------------------------------------------------------
72
73      real (kind=dbl_kind), dimension(:), allocatable, save ::
74     &     coslat, sinlat, ! cosine, sine of grid lats (for distance)
75     &     coslon, sinlon, ! cosine, sine of grid lons (for distance)
76     &     wgtstmp         ! an array to hold the link weight
77
78!***********************************************************************
79
80      contains
81
82!***********************************************************************
83
84      subroutine remap_gauswgt (lextrapdone, num_neighbors, rl_varmul)
85
86!-----------------------------------------------------------------------
87!
88!     this routine computes the inverse-distance weights for a
89!     nearest-neighbor interpolation.
90!
91!-----------------------------------------------------------------------
92
93      LOGICAL ::
94     &           lextrapdone   ! logical, true if EXTRAP done on field
95
96      REAL (kind=real_kind) ::
97     $    rl_varmul             ! Gaussian variance
98
99      INTEGER (kind=int_kind) ::
100     &       num_neighbors     ! number of neighbours
101
102!-----------------------------------------------------------------------
103!
104!     local variables
105!
106!-----------------------------------------------------------------------
107
108      logical (kind=log_kind), dimension(num_neighbors) ::
109     &     nbr_mask        ! mask at nearest neighbors
110
111      integer (kind=int_kind) :: n,
112     &     dst_add,        ! destination address
113     &     nmap            ! index of current map being computed
114
115      integer (kind=int_kind), dimension(num_neighbors) ::
116     &     nbr_add         ! source address at nearest neighbors
117
118      real (kind=dbl_kind), dimension(num_neighbors) ::
119     &     nbr_dist        ! angular distance four nearest neighbors
120
121      real (kind=dbl_kind) ::
122     &     coslat_dst,     ! cos(lat) of destination grid point
123     &     coslon_dst,     ! cos(lon) of destination grid point
124     &     sinlat_dst,     ! sin(lat) of destination grid point
125     &     sinlon_dst,     ! sin(lon) of destination grid point
126     &     dist_tot,       ! sum of neighbor distances (for normalizing)
127     &     dist_average    ! average distance between the source points
128
129!-----------------------------------------------------------------------
130!
131!     compute mappings from grid1 to grid2
132!
133!-----------------------------------------------------------------------
134
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      !*** compute the average of the distances between the source points
158      !***
159
160      call grid_dist_average(grid1_size,
161     &                       coslat, coslon,
162     &                       sinlat, sinlon,
163     &                       dist_average)
164      !***
165      !*** loop over destination grid
166      !***
167
168      grid_loop1: do dst_add = 1, grid2_size
169
170        if (.not. grid2_mask(dst_add)) cycle grid_loop1
171
172        coslat_dst = cos(grid2_center_lat(dst_add))
173        coslon_dst = cos(grid2_center_lon(dst_add))
174        sinlat_dst = sin(grid2_center_lat(dst_add))
175        sinlon_dst = sin(grid2_center_lon(dst_add))
176
177        !***
178        !*** find nearest grid points on source grid and
179        !*** distances to each point
180        !***
181
182        call grid_search_nbrg(nbr_add, nbr_dist, 
183     &                        grid2_center_lat(dst_add),
184     &                        grid2_center_lon(dst_add),
185     &                        coslat_dst, coslon_dst, 
186     &                        sinlat_dst, sinlon_dst,
187     &                        bin_addr1, 
188     &                        dist_average,  num_neighbors, rl_varmul)
189
190        !***
191        !*** compute weights based on inverse distance
192        !*** if mask is false, eliminate those points
193        !***
194
195        dist_tot = zero
196        do n=1,num_neighbors
197          if ((grid1_mask(nbr_add(n))) .or.
198     &        (.not. grid1_mask(nbr_add(n)) .and. lextrapdone)) THEN
199            nbr_dist(n) = one/nbr_dist(n)
200            dist_tot = dist_tot + nbr_dist(n)
201            nbr_mask(n) = .true.
202          else
203            nbr_mask(n) = .false.
204          endif
205        end do
206
207        !***
208        !*** normalize weights and store the link
209        !***
210
211        do n=1,num_neighbors
212          if (nbr_mask(n)) then
213            wgtstmp(1) = nbr_dist(n)/dist_tot
214            call store_link_nbr(nbr_add(n), dst_add, wgtstmp, nmap)
215            grid2_frac(dst_add) = one
216          endif
217        end do
218
219      end do grid_loop1
220
221      deallocate (coslat, coslon, sinlat, sinlon)
222
223!-----------------------------------------------------------------------
224!
225!     compute mappings from grid2 to grid1 if necessary
226!
227!-----------------------------------------------------------------------
228
229      if (num_maps > 1) then
230
231      nmap = 2
232
233      !***
234      !*** compute cos, sin of lat/lon on source grid for distance
235      !*** calculations
236      !***
237
238      allocate (coslat(grid2_size), coslon(grid2_size),
239     &          sinlat(grid2_size), sinlon(grid2_size))
240
241      coslat = cos(grid2_center_lat)
242      coslon = cos(grid2_center_lon)
243      sinlat = sin(grid2_center_lat)
244      sinlon = sin(grid2_center_lon)
245
246      !***
247      !*** compute the average of the distances between the source points
248      !***
249
250      call grid_dist_average(grid2_size,
251     &                       coslat, coslon,
252     &                       sinlat, sinlon,
253     &                       dist_average)
254
255      !***
256      !*** loop over destination grid
257      !***
258
259      grid_loop2: do dst_add = 1, grid1_size
260
261        if (.not. grid1_mask(dst_add)) cycle grid_loop2
262
263        coslat_dst = cos(grid1_center_lat(dst_add))
264        coslon_dst = cos(grid1_center_lon(dst_add))
265        sinlat_dst = sin(grid1_center_lat(dst_add))
266        sinlon_dst = sin(grid1_center_lon(dst_add))
267
268        !***
269        !*** find four nearest grid points on source grid and
270        !*** distances to each point
271        !***
272
273        call grid_search_nbrg(nbr_add, nbr_dist,
274     &                        grid1_center_lat(dst_add),
275     &                        grid1_center_lon(dst_add),
276     &                        coslat_dst, coslon_dst, 
277     &                        sinlat_dst, sinlon_dst,
278     &                        bin_addr2, 
279     &                        dist_average, num_neighbors, rl_varmul)
280
281        !***
282        !*** compute weights based on inverse distance
283        !*** if mask is false, eliminate those points
284        !***
285
286        dist_tot = zero
287        do n=1,num_neighbors
288          if (grid2_mask(nbr_add(n))) then
289            nbr_dist(n) = one/nbr_dist(n)
290            dist_tot = dist_tot + nbr_dist(n)
291            nbr_mask(n) = .true.
292          else
293            nbr_mask(n) = .false.
294          endif
295        end do
296
297        !***
298        !*** normalize weights and store the link
299        !***
300
301        do n=1,num_neighbors
302          if (nbr_mask(n)) then
303            wgtstmp(1) = nbr_dist(n)/dist_tot
304            call store_link_nbr(dst_add, nbr_add(n), wgtstmp, nmap)
305            grid1_frac(dst_add) = one
306          endif
307        end do
308
309      end do grid_loop2
310
311      deallocate (coslat, coslon, sinlat, sinlon)
312
313      endif
314
315      deallocate(wgtstmp)
316
317!-----------------------------------------------------------------------
318
319      end subroutine remap_gauswgt
320
321!***********************************************************************
322
323      subroutine grid_search_nbrg(nbr_add, nbr_dist, plat, plon, 
324     &               coslat_dst, coslon_dst, sinlat_dst, sinlon_dst,
325     &               src_bin_add, dst_average,
326     $               num_neighbors, rl_varmul)
327
328!-----------------------------------------------------------------------
329!
330!     this routine finds the closest num_neighbor points to a search
331!     point and computes a distance to each of the neighbors.
332!
333!-----------------------------------------------------------------------
334
335!-----------------------------------------------------------------------
336!
337!     input variables
338!
339!-----------------------------------------------------------------------
340
341      integer (kind=int_kind), dimension(:,:), intent(in) ::
342     &        src_bin_add ! search bins for restricting search
343
344      real (kind=dbl_kind), intent(in) ::
345     &        plat,         ! latitude  of the search point
346     &        plon,         ! longitude of the search point
347     &        coslat_dst,   ! cos(lat)  of the search point
348     &        coslon_dst,   ! cos(lon)  of the search point
349     &        sinlat_dst,   ! sin(lat)  of the search point
350     &        sinlon_dst,   ! sin(lon)  of the search point
351     &        dst_average   ! average distance between the source points
352
353      REAL (kind=real_kind), intent(in) ::
354     &        rl_varmul     ! Gaussian variance
355
356      INTEGER (kind=int_kind) ::
357     &       num_neighbors     ! number of neighbours
358
359!-----------------------------------------------------------------------
360!
361!     output variables
362!
363!-----------------------------------------------------------------------
364
365      integer (kind=int_kind), dimension(num_neighbors), intent(out) ::
366     &        nbr_add  ! address of each of the closest points
367
368      real (kind=dbl_kind), dimension(num_neighbors), intent(out) ::
369     &        nbr_dist ! distance to each of the closest points
370
371!-----------------------------------------------------------------------
372!
373!     local variables
374!
375!-----------------------------------------------------------------------
376
377      integer (kind=int_kind) :: n, nmax, nadd, nchk, ! dummy indices
378     &        min_add, max_add, nm1, np1, i, j, ip1, im1, jp1, jm1
379
380      real (kind=dbl_kind) ::
381     &        distance      ! angular distance
382
383      real (kind=dbl_kind) ::
384     &        variance      ! variance for the gaussian FUNCTION
385
386
387!-----------------------------------------------------------------------
388!
389!     loop over source grid and find nearest neighbors
390!
391!-----------------------------------------------------------------------
392
393      !***
394      !*** restrict the search using search bins
395      !*** expand the bins to catch neighbors
396      !***
397
398      select case (restrict_type)
399      case('LATITUDE')
400
401        do n=1,num_srch_bins
402          if (plat >= bin_lats(1,n) .and. plat <= bin_lats(2,n)) then
403            min_add = src_bin_add(1,n)
404            max_add = src_bin_add(2,n)
405
406            nm1 = max(n-1,1)
407            np1 = min(n+1,num_srch_bins)
408
409            min_add = min(min_add,src_bin_add(1,nm1))
410            max_add = max(max_add,src_bin_add(2,nm1))
411            min_add = min(min_add,src_bin_add(1,np1))
412            max_add = max(max_add,src_bin_add(2,np1))
413          endif
414        end do
415
416      case('LATLON')
417
418        n = 0
419        nmax = nint(sqrt(real(num_srch_bins)))
420        do j=1,nmax
421        jp1 = min(j+1,nmax)
422        jm1 = max(j-1,1)
423        do i=1,nmax
424          ip1 = min(i+1,nmax)
425          im1 = max(i-1,1)
426
427          n = n+1
428          if (plat >= bin_lats(1,n) .and. plat <= bin_lats(2,n) .and.
429     &        plon >= bin_lons(1,n) .and. plon <= bin_lons(2,n)) then
430            min_add = src_bin_add(1,n)
431            max_add = src_bin_add(2,n)
432
433            nm1 = (jm1-1)*nmax + im1
434            np1 = (jp1-1)*nmax + ip1
435            nm1 = max(nm1,1)
436            np1 = min(np1,num_srch_bins)
437
438            min_add = min(min_add,src_bin_add(1,nm1))
439            max_add = max(max_add,src_bin_add(2,nm1))
440            min_add = min(min_add,src_bin_add(1,np1))
441            max_add = max(max_add,src_bin_add(2,np1))
442          endif
443        end do
444        end do
445
446      end select
447
448      !***
449      !*** initialize distance and address arrays
450      !***
451
452      nbr_add = 0
453      nbr_dist = bignum
454      variance = rl_varmul*dst_average*dst_average
455      do nadd=min_add,max_add
456
457        !***
458        !*** find distance to this point
459        !***
460
461        distance = acos(sinlat_dst*sinlat(nadd) +
462     &                  coslat_dst*coslat(nadd)*
463     &                 (coslon_dst*coslon(nadd) +
464     &                  sinlon_dst*sinlon(nadd)) )
465        !distance = min(500., distance) !ts-sv
466        !distance = exp(.5*distance*distance/variance)
467        !***
468        !*** store the address and distance if this is one of the
469        !*** smallest four so far
470        !***
471
472        check_loop: do nchk=1,num_neighbors
473          if (distance .lt. nbr_dist(nchk)) THEN
474            do n=num_neighbors,nchk+1,-1
475              nbr_add(n) = nbr_add(n-1)
476              nbr_dist(n) = nbr_dist(n-1)
477            end do
478            nbr_add(nchk) = nadd
479            nbr_dist(nchk) = distance
480            exit check_loop
481          endif
482        end do check_loop
483
484      end do
485
486      exp_loop: do nchk=1,num_neighbors
487          nbr_dist(nchk) = 
488     &     exp(.5*nbr_dist(nchk)*nbr_dist(nchk)/variance)
489      end do exp_loop
490
491!-----------------------------------------------------------------------
492
493      end subroutine grid_search_nbrg 
494
495!***********************************************************************
496
497      subroutine store_link_nbr(add1, add2, weights, nmap)
498
499!-----------------------------------------------------------------------
500!
501!     this routine stores the address and weight for this link in
502!     the appropriate address and weight arrays and resizes those
503!     arrays if necessary.
504!
505!-----------------------------------------------------------------------
506
507!-----------------------------------------------------------------------
508!
509!     input variables
510!
511!-----------------------------------------------------------------------
512
513      integer (kind=int_kind), intent(in) ::
514     &        add1,  ! address on grid1
515     &        add2,  ! address on grid2
516     &        nmap   ! identifies which direction for mapping
517
518      real (kind=dbl_kind), dimension(:), intent(in) ::
519     &        weights ! array of remapping weights for this link
520
521!-----------------------------------------------------------------------
522!
523!     increment number of links and check to see if remap arrays need
524!     to be increased to accomodate the new link.  then store the
525!     link.
526!
527!-----------------------------------------------------------------------
528
529      select case (nmap)
530      case(1)
531
532        num_links_map1  = num_links_map1 + 1
533
534        if (num_links_map1 > max_links_map1) 
535     &     call resize_remap_vars(1,resize_increment)
536
537        grid1_add_map1(num_links_map1) = add1
538        grid2_add_map1(num_links_map1) = add2
539        wts_map1    (:,num_links_map1) = weights
540
541      case(2)
542
543        num_links_map2  = num_links_map2 + 1
544
545        if (num_links_map2 > max_links_map2) 
546     &     call resize_remap_vars(2,resize_increment)
547
548        grid1_add_map2(num_links_map2) = add1
549        grid2_add_map2(num_links_map2) = add2
550        wts_map2    (:,num_links_map2) = weights
551
552      end select
553
554!-----------------------------------------------------------------------
555
556      end subroutine store_link_nbr
557
558!***********************************************************************
559
560      subroutine grid_dist_average(grid_size,
561     &                             coslat_grid, coslon_grid,
562     &                             sinlat_grid, sinlon_grid,
563     &                             dist_average)
564
565!-----------------------------------------------------------------------
566!
567!     this routine computes the average distance between the points of a
568!     grid.
569!
570!-----------------------------------------------------------------------
571
572!-----------------------------------------------------------------------
573!
574!     output variables
575!
576!-----------------------------------------------------------------------
577
578      real (kind=dbl_kind), intent(out) ::
579     &        dist_average ! distance to each of the closest points
580
581!-----------------------------------------------------------------------
582!
583!     input variables
584!
585!-----------------------------------------------------------------------
586
587      integer (kind=int_kind), intent(in) ::
588     &        grid_size 
589
590      real (kind=dbl_kind), dimension(:), intent(in) ::
591     &        coslat_grid,   ! cos(lat)  of the grid points
592     &        coslon_grid,   ! cos(lon)  of the grid points
593     &        sinlat_grid,   ! sin(lat)  of the grid points
594     &        sinlon_grid    ! sin(lon)  of the grid points
595
596!-----------------------------------------------------------------------
597!
598!     local variables
599!
600!-----------------------------------------------------------------------
601
602      integer (kind=int_kind) :: i
603
604!-----------------------------------------------------------------------
605!
606!     compute the distance over the grid and average
607!
608!-----------------------------------------------------------------------
609
610      dist_average = 0.0
611      DO i = 1, grid_size
612        IF (i .eq. 1) THEN
613            dist_average = dist_average +
614     &          acos(sinlat_grid(grid_size)*sinlat_grid(i) +
615     &          coslat_grid(grid_size)*coslat_grid(i)*
616     &          (coslon_grid(grid_size)*coslon_grid(i) +
617     &          sinlon_grid(grid_size)*sinlon_grid(i))) 
618            dist_average = dist_average +
619     &          acos(sinlat_grid(i)*sinlat_grid(i+1) +
620     &          coslat_grid(i)*coslat_grid(i+1)*
621     &          (coslon_grid(i)*coslon_grid(i+1) +
622     &          sinlon_grid(i)*sinlon_grid(i+1)))             
623        ELSE IF (i .eq. grid_size) THEN
624            dist_average = dist_average +
625     &          acos(sinlat_grid(i-1)*sinlat_grid(i) +
626     &          coslat_grid(i-1)*coslat_grid(i)*
627     &          (coslon_grid(i-1)*coslon_grid(i) +
628     &          sinlon_grid(i-1)*sinlon_grid(i)))             
629            dist_average = dist_average +
630     &          acos(sinlat_grid(i)*sinlat_grid(1) +
631     &          coslat_grid(i)*coslat_grid(1)*
632     &          (coslon_grid(i)*coslon_grid(1) +
633     &          sinlon_grid(i)*sinlon_grid(1)))             
634        ELSE
635            dist_average = dist_average +
636     &          acos(sinlat_grid(i-1)*sinlat_grid(i) +
637     &          coslat_grid(i-1)*coslat_grid(i)*
638     &          (coslon_grid(i-1)*coslon_grid(i) +
639     &          sinlon_grid(i-1)*sinlon_grid(i)))             
640            dist_average = dist_average +
641     &          acos(sinlat_grid(i)*sinlat_grid(i+1) +
642     &          coslat_grid(i)*coslat_grid(i+1)*
643     &          (coslon_grid(i)*coslon_grid(i+1) +
644     &          sinlon_grid(i)*sinlon_grid(i+1)))             
645        END IF
646      END DO
647      dist_average = dist_average / (2*grid_size)
648
649      END subroutine grid_dist_average
650
651!***********************************************************************
652
653      end module remap_gaussian_weight
654
655!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Note: See TracBrowser for help on using the repository browser.