source: CPL/oasis3/trunk/src/lib/scrip/src/remap_bilinear_reduced.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: 22.1 KB
Line 
1C****
2C               *****************************
3C               * OASIS MODULE  -  LEVEL ? *
4C               * -------------     ------- *
5C               *****************************
6C
7C**** remap_bilinear_reduced - calculate reduced grid bilinear remapping
8C
9C     Purpose:
10C     -------
11C     Adaptation of SCRIP 1.4 remap_bilinear module to calculate 
12C     bilinear remapping for reduced grids.
13C
14C**   Interface:
15C     ---------
16C       *CALL*  *remap_bilin_reduced*
17C
18C     Input:
19C     -----
20C
21C     Output:
22C     ------
23C
24C     History:
25C     -------
26C       Version   Programmer     Date        Description
27C       -------   ----------     ----        ----------- 
28C       2.5       D. Declat      2002/07     created
29C       2.5       S. Valcke      2002/09     modified
30C
31C %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
32!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
33!
34!     this module contains necessary routines for performing an
35!     bilinear interpolation on reduced grids.
36!
37!-----------------------------------------------------------------------
38!
39!     CVS:$Id: remap_bilinear_reduced.f,v 1.1.1.1 2005/03/23 16:01:11 adm Exp $
40!
41!     Copyright (c) 1997, 1998 the Regents of the University of
42!       California.
43!
44!     This software and ancillary information (herein called software)
45!     called SCRIP is made available under the terms described here. 
46!     The software has been approved for release with associated
47!     LA-CC Number 98-45.
48!
49!     Unless otherwise indicated, this software has been authored
50!     by an employee or employees of the University of California,
51!     operator of the Los Alamos National Laboratory under Contract
52!     No. W-7405-ENG-36 with the U.S. Department of Energy.  The U.S.
53!     Government has rights to use, reproduce, and distribute this
54!     software.  The public may copy and use this software without
55!     charge, provided that this Notice and any statement of authorship
56!     are reproduced on all copies.  Neither the Government nor the
57!     University makes any warranty, express or implied, or assumes
58!     any liability or responsibility for the use of this software.
59!
60!     If software is modified to produce derivative works, such modified
61!     software should be clearly marked, so as not to confuse it with
62!     the version available from Los Alamos National Laboratory.
63!
64!***********************************************************************
65
66      module remap_bilinear_reduced
67
68!-----------------------------------------------------------------------
69
70      use kinds_mod     ! defines common data types
71      use constants     ! defines common constants
72      use grids         ! module containing grid info
73      use remap_vars    ! module containing remap info
74      USE mod_unit
75      USE mod_printing
76
77      implicit none
78
79!-----------------------------------------------------------------------
80
81      integer (kind=int_kind), parameter ::
82     &    max_iter = 100   ! max iteration count for i,j iteration
83
84      real (kind=dbl_kind), parameter ::
85     &     converge = 1.e-10_dbl_kind  ! convergence criterion
86
87!***********************************************************************
88
89      contains
90
91!***********************************************************************
92
93      subroutine remap_bilin_reduced (lextrapdone)
94
95!-----------------------------------------------------------------------
96!
97!     this routine computes the weights for a bilinear interpolation.
98!
99!-----------------------------------------------------------------------
100
101      LOGICAL ::
102     &           lextrapdone   ! logical, true if EXTRAP done on field
103
104!-----------------------------------------------------------------------
105!
106!     local variables
107!
108!-----------------------------------------------------------------------
109
110      integer (kind=int_kind) :: n, icount, i,
111     &     dst_add,        ! destination address
112     &     iter,           ! iteration counter
113     &     nmap           ! index of current map being computed
114
115      integer (kind=int_kind), dimension(4) :: 
116     &     src_add         ! address for the four source points
117
118      real (kind=dbl_kind), dimension(4)  ::
119     &     src_lats,       ! latitudes  of four bilinear corners
120     &     src_lons,       ! longitudes of four bilinear corners
121     &     wgts            ! bilinear weights for four corners
122
123      real (kind=dbl_kind) ::
124     &     plat, plon,       ! lat/lon coords of destination point
125     &     iguess, jguess,   ! current guess for bilinear coordinate
126     &     deli, delj,       ! corrections to i,j
127     &     dth1, dth2, dth3, ! some latitude  differences
128     &     dph1, dph2, dph3, ! some longitude differences
129     &     dthp, dphp,       ! difference between point and sw corner
130     &     mat1, mat2, mat3, mat4, ! matrix elements
131     &     determinant,       ! matrix determinant
132     &     sum_wgts          ! sum of weights for normalization
133
134      real (kind=dbl_kind) :: 
135     &      coslat_dst, sinlat_dst, coslon_dst, sinlon_dst,
136     &      dist_min, distance, ! for computing dist-weighted avg
137     &      src_latsnn
138
139      integer (kind=int_kind) :: min_add, max_add, srch_add, src_addnn
140
141!-----------------------------------------------------------------------
142!
143!     compute mappings from grid1 to grid2
144!
145!-----------------------------------------------------------------------
146
147      nmap = 1
148      if (grid1_rank /= 2) then
149          stop 'Can not do bilinear interpolation when grid_rank /= 2'
150      endif
151
152      !***
153      !*** loop over destination grid
154      !***
155
156      grid_loop1: do dst_add = 1, grid2_size
157
158         if (.not. grid2_mask(dst_add)) cycle grid_loop1
159
160         plat = grid2_center_lat(dst_add)
161         plon = grid2_center_lon(dst_add)
162
163        !***
164        !*** find nearest square of grid points on source grid
165        !***
166
167        call grid_search_bilin_rd(src_add, src_lats, src_lons, 
168     &                            min_add, max_add,
169     &                            plat, plon, grid1_dims,
170     &                            grid1_center_lat, grid1_center_lon,
171     &                            lextrapdone)
172
173        if (src_add(1) > 0) THEN
174
175        !***
176        !*** if the 4 surrounding points have been found and are
177        !*** non-masked, do bilinear interpolation
178        !***
179
180          grid2_frac(dst_add) = one
181
182          !***
183          !*** iterate to find i,j for bilinear approximation
184          !***
185
186          dth1 = src_lats(2) - src_lats(1)
187          dth2 = src_lats(4) - src_lats(1)
188          dth3 = src_lats(3) - src_lats(2) - dth2
189
190          dph1 = src_lons(2) - src_lons(1)
191          dph2 = src_lons(4) - src_lons(1)
192          dph3 = src_lons(3) - src_lons(2)
193
194          if (dph1 >  three*pih) dph1 = dph1 - pi2
195          if (dph2 >  three*pih) dph2 = dph2 - pi2
196          if (dph3 >  three*pih) dph3 = dph3 - pi2
197          if (dph1 < -three*pih) dph1 = dph1 + pi2
198          if (dph2 < -three*pih) dph2 = dph2 + pi2
199          if (dph3 < -three*pih) dph3 = dph3 + pi2
200
201          dph3 = dph3 - dph2
202
203          iguess = half
204          jguess = half
205
206          iter_loop1: do iter=1,max_iter
207
208            dthp = plat - src_lats(1) - dth1*iguess -
209     &                    dth2*jguess - dth3*iguess*jguess
210            dphp = plon - src_lons(1)
211
212            if (dphp >  three*pih) dphp = dphp - pi2
213            if (dphp < -three*pih) dphp = dphp + pi2
214
215            dphp = dphp - dph1*iguess - dph2*jguess - 
216     &                    dph3*iguess*jguess
217
218            mat1 = dth1 + dth3*jguess
219            mat2 = dth2 + dth3*iguess
220            mat3 = dph1 + dph3*jguess
221            mat4 = dph2 + dph3*iguess
222
223            determinant = mat1*mat4 - mat2*mat3
224
225            deli = (dthp*mat4 - mat2*dphp)/determinant
226            delj = (mat1*dphp - dthp*mat3)/determinant
227
228            if (abs(deli) < converge .and. 
229     &          abs(delj) < converge) exit iter_loop1
230
231            iguess = iguess + deli
232            jguess = jguess + delj
233
234          end do iter_loop1
235
236          if (iter <= max_iter) then
237
238            !***
239            !*** successfully found i,j - compute weights
240            !***
241
242            wgts(1) = (one-iguess)*(one-jguess)
243            wgts(2) = iguess*(one-jguess)
244            wgts(3) = iguess*jguess
245            wgts(4) = (one-iguess)*jguess
246            call store_link_bilin(dst_add, src_add, wgts, nmap)
247          else
248            print *,'Point coords: ',plat,plon
249            print *,'Dest grid lats: ',src_lats
250            print *,'Dest grid lons: ',src_lons
251            print *,'Dest grid addresses: ',src_add
252            print *,'Current i,j : ',iguess, jguess
253            stop 'Iteration for i,j exceed max iteration count'
254          endif
255
256        else if (src_add(1) < 0) THEN
257
258          !***
259          !*** We are in the first or last bin or at least one of the 4
260          !*** neighbours was masked. Do distance-weighted average using
261          !*** the non-masked points among the 4 closest ones.
262
263          IF (nlogprt .eq. 2) then
264              WRITE(nulou,*) ' '
265              WRITE(nulou,*) 
266     &   'WARNING: Cannot make bilinear interpolation for target point'
267     &            ,dst_add
268              WRITE(nulou,*) 
269     &    'Using non-masked points among 4 nearest neighbors.'
270              WRITE(nulou,*) ' '
271          ENDIF
272
273          !***
274          !*** Find the 4 closest points
275          !***
276
277          coslat_dst = cos(plat)
278          sinlat_dst = sin(plat)
279          coslon_dst = cos(plon)
280          sinlon_dst = sin(plon)
281          src_add = 0
282          dist_min = bignum
283          src_lats = bignum
284          do srch_add = min_add,max_add
285            distance=acos(coslat_dst*cos(grid1_center_lat(srch_add))*
286     &          (coslon_dst*cos(grid1_center_lon(srch_add)) +
287     &          sinlon_dst*sin(grid1_center_lon(srch_add)))+
288     &          sinlat_dst*sin(grid1_center_lat(srch_add)))
289
290            if (distance < dist_min) then
291                sort_loop: do n=1,4
292                if (distance < src_lats(n)) then
293                    do i=4,n+1,-1
294                      src_add (i) = src_add (i-1)
295                      src_lats(i) = src_lats(i-1)
296                    end do
297                    src_add (n) = srch_add
298                    src_lats(n) = distance
299                    dist_min = src_lats(4)
300                    exit sort_loop
301                endif
302                end do sort_loop
303            endif
304          end do
305
306          src_lons = one/(src_lats + tiny)
307          distance = sum(src_lons)
308          src_lats = src_lons/distance
309
310          !***
311          !*** Among 4 closest points, keep only the non-masked ones
312          !***
313
314          icount = 0
315          do n=1,4
316            if (grid1_mask(src_add(n)) .or.
317     &          (.not. grid1_mask(src_add(n)) .and. lextrapdone)) then
318                icount = icount + 1
319            else
320                src_lats(n) = zero
321            endif
322          end do
323
324          if (icount > 0) then
325              !*** renormalize weights
326              sum_wgts = sum(src_lats)
327              wgts(1) = src_lats(1)/sum_wgts
328              wgts(2) = src_lats(2)/sum_wgts
329              wgts(3) = src_lats(3)/sum_wgts
330              wgts(4) = src_lats(4)/sum_wgts
331
332              grid2_frac(dst_add) = one
333              call store_link_bilin(dst_add, src_add, wgts, nmap)
334          ELSE
335              WRITE(nulou,*) '  '
336              WRITE(nulou,*) 
337     &       'All 4 surrounding source grid points are masked'
338              WRITE(nulou,*) 'for target point ',dst_add
339              WRITE(nulou,*) 'with longitude and latitude', plon, plat
340              WRITE(nulou,*) 'Using the nearest non-masked neighbour.' 
341              WRITE(nulou,*) '  '
342              src_latsnn = bignum
343              do srch_add = min_add,max_add
344                if (grid1_mask(srch_add) .or.
345     &          (.not. grid1_mask(srch_add) .and. lextrapdone)) then
346                    distance=acos(coslat_dst
347     &                  *cos(grid1_center_lat(srch_add))*
348     &                  (coslon_dst*cos(grid1_center_lon(srch_add)) +
349     &                  sinlon_dst*sin(grid1_center_lon(srch_add)))+
350     &                  sinlat_dst*sin(grid1_center_lat(srch_add)))
351
352                    if (distance < src_latsnn) then
353                        src_addnn = srch_add
354                        src_latsnn = distance
355                    endif
356                endif
357              end DO
358
359              WRITE(nulou,*) 
360     &        'Nearest non masked neighbour is source point ',src_addnn
361              WRITE(nulou,*) 'with longitude and latitude',
362     &        grid1_center_lon(src_addnn), grid1_center_lat(src_addnn) 
363              WRITE(nulou,*) '  '
364
365              wgts(1) = 1.
366              wgts(2) = 0.
367              wgts(3) = 0.
368              wgts(4) = 0.
369              src_add(1) = src_addnn
370              src_add(2) = 0
371              src_add(3) = 0
372              src_add(4) = 0
373
374              grid2_frac(dst_add) = one
375              call store_link_bilin(dst_add, src_add, wgts, nmap)
376
377          ENDIF
378      ENDIF
379      end do grid_loop1
380!
381!-----------------------------------------------------------------------
382
383      end subroutine remap_bilin_reduced
384
385!***********************************************************************
386
387      subroutine grid_search_bilin_rd(src_add, src_lats, src_lons,
388     &                                min_add, max_add,
389     &                                plat, plon, src_grid_dims,
390     &                                src_center_lat, src_center_lon,
391     &                                lextrapdone)
392
393!-----------------------------------------------------------------------
394!
395!     this routine finds the location of the search point plat, plon
396!     in the source grid and returns the corners needed for a bilinear
397!     interpolation. The target grid is a reduced grid.
398!
399!-----------------------------------------------------------------------
400
401!-----------------------------------------------------------------------
402!
403!     output variables
404!
405!-----------------------------------------------------------------------
406
407      integer (kind=int_kind), dimension(4), intent(out) ::
408     &        src_add  ! address of each corner point enclosing P
409
410      real (kind=dbl_kind), dimension(4), intent(out) ::
411     &        src_lats, ! latitudes  of the four corner points
412     &        src_lons  ! longitudes of the four corner points
413
414      integer (kind=int_kind) :: min_add, max_add
415
416!-----------------------------------------------------------------------
417!
418!     input variables
419!
420!-----------------------------------------------------------------------
421
422      real (kind=dbl_kind), intent(in) ::
423     &        plat,   ! latitude  of the search point
424     &        plon    ! longitude of the search point
425
426      integer (kind=int_kind), dimension(2), intent(in) ::
427     &        src_grid_dims  ! size of each src grid dimension
428
429      real (kind=dbl_kind), dimension(:), intent(in) ::
430     &        src_center_lat, ! latitude  of each src grid center
431     &        src_center_lon  ! longitude of each src grid center
432
433      LOGICAL ::
434     &           lextrapdone   ! logical, true if EXTRAP done on field
435
436!-----------------------------------------------------------------------
437!
438!     local variables
439!
440!-----------------------------------------------------------------------
441
442      integer (kind=int_kind) :: n, next_n, srch_add, ni,  ! dummy indices
443     &    nx, ny, ntotmask,           ! dimensions of src grid
444     &    inter_add ! add for restricting search
445!
446      integer (kind=int_kind), DIMENSION(4) :: src_bid
447
448!-----------------------------------------------------------------------
449!
450!     restrict search first using bins
451!
452!-----------------------------------------------------------------------
453
454      nx = src_grid_dims(1)
455      inter_add = 0
456
457      src_add = 0
458
459      min_add = size(src_center_lat) + 1
460      max_add = 1
461      if (plat >= bin_lats_r(1,1)) then
462          min_add = 0
463          max_add = bin_addr1_r(4,1)
464          inter_add = bin_addr1_r(3,1)
465      else if (plat <= bin_lats_r(1,num_srch_red)) then
466          max_add = nx + 1
467          min_add = bin_addr1_r(1,num_srch_red)
468          inter_add = bin_addr1_r(3,num_srch_red)
469      else
470          srch_loopred: do n=1,num_srch_red
471            if (plat <= bin_lats_r(1,n) 
472     &        .and. plat >= bin_lats_r(2,n)) then
473                min_add = bin_addr1_r(1,n)
474                max_add = bin_addr1_r(4,n)
475                inter_add = bin_addr1_r(3,n)
476                exit srch_loopred
477            endif
478          end DO srch_loopred
479      ENDIF
480
481!-----------------------------------------------------------------------
482!
483!     now perform a more detailed search
484!
485!-----------------------------------------------------------------------
486
487      if (min_add .ne. 0 .and. max_add .ne. nx+1) THEN
488         !*** The concerned bins are not the top north or south ones.
489         !*** We should be able to find the four corners
490         !*** for the bilinear interpolation.
491
492          IF ( plon <= src_center_lon(min_add) ) THEN
493              src_add(1) = inter_add-1
494              src_add(2) = min_add
495          ELSE IF ( plon > src_center_lon(inter_add-1) ) then
496              src_add(1) = inter_add-1
497              src_add(2) = min_add
498          else
499              srch_loop2: do srch_add = min_add, inter_add-2
500                 if ( (plon > src_center_lon(srch_add)) .and.
501     &            (plon <= src_center_lon(srch_add+1)) )then
502                     src_add(1) = srch_add
503                     src_add(2) = srch_add+1
504                     exit srch_loop2
505                 endif
506               end do srch_loop2
507           ENDIF
508           IF ( plon <= src_center_lon(inter_add) ) THEN
509               src_add(3) = max_add
510               src_add(4) = inter_add
511           ELSE IF ( plon >= src_center_lon(max_add) ) then
512               src_add(3) = max_add
513               src_add(4) = inter_add
514           else
515               srch_loop3: do srch_add = inter_add, max_add
516                 if ( (plon >= src_center_lon(srch_add)) .and.
517     &             (plon <= src_center_lon(srch_add+1)) )then
518                      src_add(3) = srch_add
519                      src_add(4) = srch_add+1
520                      exit srch_loop3
521                  endif
522               enddo srch_loop3
523           ENDIF
524           src_lats(1) = src_center_lat(src_add(3))
525           src_lats(2) = src_center_lat(src_add(4))
526           src_lats(3) = src_center_lat(src_add(2))
527           src_lats(4) = src_center_lat(src_add(1))
528     
529           src_lons(1) = src_center_lon(src_add(3))
530           src_lons(2) = src_center_lon(src_add(4))
531           src_lons(3) = src_center_lon(src_add(2))
532           src_lons(4) = src_center_lon(src_add(1))
533
534           src_bid=src_add
535
536           src_add(1) = src_bid(3)
537           src_add(2) = src_bid(4)
538           src_add(3) = src_bid(2)
539           src_add(4) = src_bid(1)
540   
541           ! Check if one point is masked; IF so, nearest-neighbour
542           ! interpolation will be used
543
544           ntotmask = 0
545           do ni=1,4
546             if (.not. grid1_mask(src_add(ni)).and. 
547     &           .not. lextrapdone) 
548     &           ntotmask = ntotmask + 1 
549           end DO
550           IF (ntotmask .gt. 0) src_add(1) = -src_add(1) 
551         
552       ELSE 
553
554           !*** We are in the first or last bin.  Put src_add = -1 so that
555           !*** distance-weighted average of the 4 non-masked closest points
556           !*** is done in calling routine.
557
558           src_add = -1
559
560       ENDIF
561
562!-----------------------------------------------------------------------
563
564      end subroutine grid_search_bilin_rd 
565
566!***********************************************************************
567
568      subroutine store_link_bilin(dst_add, src_add, weights, nmap)
569
570!-----------------------------------------------------------------------
571!
572!     this routine stores the address and weight for four links
573!     associated with one destination point in the appropriate address
574!     and weight arrays and resizes those arrays if necessary.
575!
576!-----------------------------------------------------------------------
577
578!-----------------------------------------------------------------------
579!
580!     input variables
581!
582!-----------------------------------------------------------------------
583
584      integer (kind=int_kind), intent(in) ::
585     &        dst_add,  ! address on destination grid
586     &        nmap      ! identifies which direction for mapping
587
588      integer (kind=int_kind), dimension(4), intent(in) ::
589     &        src_add   ! addresses on source grid
590
591      real (kind=dbl_kind), dimension(4), intent(in) ::
592     &        weights ! array of remapping weights for these links
593
594!-----------------------------------------------------------------------
595!
596!     local variables
597!
598!-----------------------------------------------------------------------
599
600      integer (kind=int_kind) :: n, ! dummy index
601     &       num_links_old          ! placeholder for old link number
602
603!-----------------------------------------------------------------------
604!
605!     increment number of links and check to see if remap arrays need
606!     to be increased to accomodate the new link.  then store the
607!     link.
608!
609!-----------------------------------------------------------------------
610
611      select case (nmap)
612      case(1)
613
614        num_links_old  = num_links_map1
615        num_links_map1 = num_links_old + 4
616
617        if (num_links_map1 > max_links_map1) 
618     &     call resize_remap_vars(1,resize_increment)
619
620        do n=1,4
621          grid1_add_map1(num_links_old+n) = src_add(n)
622          grid2_add_map1(num_links_old+n) = dst_add
623          wts_map1    (1,num_links_old+n) = weights(n)
624        end do
625
626      case(2)
627
628        num_links_old  = num_links_map2
629        num_links_map2 = num_links_old + 4
630
631        if (num_links_map2 > max_links_map2) 
632     &     call resize_remap_vars(2,resize_increment)
633
634        do n=1,4
635          grid1_add_map2(num_links_old+n) = dst_add
636          grid2_add_map2(num_links_old+n) = src_add(n)
637          wts_map2    (1,num_links_old+n) = weights(n)
638        end do
639
640      end select
641
642!-----------------------------------------------------------------------
643
644      end subroutine store_link_bilin
645
646!***********************************************************************
647
648      end module remap_bilinear_reduced
649
650!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Note: See TracBrowser for help on using the repository browser.