source: CPL/oasis3-mct/branches/OASIS3-MCT_2.0_branch/lib/scrip/src/remap_bilinear_reduced.f @ 4775

Last change on this file since 4775 was 4775, checked in by aclsce, 5 years ago
  • Imported oasis3-mct from Cerfacs svn server (not suppotred anymore).

The version has been extracted from ​https://oasis3mct.cerfacs.fr/svn/branches/OASIS3-MCT_2.0_branch/oasis3-mct@1818
File size: 23.3 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 2826 2010-12-10 11:14:21Z valcke $
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_oasis_flush
75
76      implicit none
77
78!-----------------------------------------------------------------------
79
80      integer (kind=int_kind), parameter ::
81     &    max_iter = 100   ! max iteration count for i,j iteration
82
83      real (kind=dbl_kind), parameter ::
84     &     converge = 1.e-10_dbl_kind  ! convergence criterion
85
86!***********************************************************************
87
88      contains
89
90!***********************************************************************
91
92      subroutine remap_bilin_reduced (lextrapdone)
93
94!-----------------------------------------------------------------------
95!
96!     this routine computes the weights for a bilinear interpolation.
97!
98!-----------------------------------------------------------------------
99
100      LOGICAL ::
101     &           lextrapdone   ! logical, true if EXTRAP done on field
102      LOGICAL :: ll_nnei        ! true (default) if extra search is done
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, arg
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      IF (nlogprt .GE. 2) THEN
148         WRITE (UNIT = nulou,FMT = *)' '
149         WRITE (UNIT = nulou,FMT = *)
150     &       'Entering routine remap_bilin_reduced'
151         CALL OASIS_FLUSH_SCRIP(nulou)
152      ENDIF
153!
154      ll_nnei = .true.
155      nmap = 1
156      if (grid1_rank /= 2) then
157          stop 'Can not do bilinear interpolation when grid_rank /= 2'
158      endif
159
160      !***
161      !*** loop over destination grid
162      !***
163
164      grid_loop1: do dst_add = 1, grid2_size
165
166         if (.not. grid2_mask(dst_add)) cycle grid_loop1
167
168         plat = grid2_center_lat(dst_add)
169         plon = grid2_center_lon(dst_add)
170
171        !***
172        !*** find nearest square of grid points on source grid
173        !***
174
175        call grid_search_bilin_rd(src_add, src_lats, src_lons, 
176     &                            min_add, max_add,
177     &                            plat, plon, grid1_dims,
178     &                            grid1_center_lat, grid1_center_lon,
179     &                            lextrapdone)
180        if (src_add(1) > 0) THEN
181
182        !***
183        !*** if the 4 surrounding points have been found and are
184        !*** non-masked, do bilinear interpolation
185        !***
186
187          grid2_frac(dst_add) = one
188
189          !***
190          !*** iterate to find i,j for bilinear approximation
191          !***
192
193          dth1 = src_lats(2) - src_lats(1)
194          dth2 = src_lats(4) - src_lats(1)
195          dth3 = src_lats(3) - src_lats(2) - dth2
196
197          dph1 = src_lons(2) - src_lons(1)
198          dph2 = src_lons(4) - src_lons(1)
199          dph3 = src_lons(3) - src_lons(2)
200
201          if (dph1 >  three*pih) dph1 = dph1 - pi2
202          if (dph2 >  three*pih) dph2 = dph2 - pi2
203          if (dph3 >  three*pih) dph3 = dph3 - pi2
204          if (dph1 < -three*pih) dph1 = dph1 + pi2
205          if (dph2 < -three*pih) dph2 = dph2 + pi2
206          if (dph3 < -three*pih) dph3 = dph3 + pi2
207
208          dph3 = dph3 - dph2
209
210          iguess = half
211          jguess = half
212
213          iter_loop1: do iter=1,max_iter
214
215            dthp = plat - src_lats(1) - dth1*iguess -
216     &                    dth2*jguess - dth3*iguess*jguess
217            dphp = plon - src_lons(1)
218
219            if (dphp >  three*pih) dphp = dphp - pi2
220            if (dphp < -three*pih) dphp = dphp + pi2
221
222            dphp = dphp - dph1*iguess - dph2*jguess - 
223     &                    dph3*iguess*jguess
224
225            mat1 = dth1 + dth3*jguess
226            mat2 = dth2 + dth3*iguess
227            mat3 = dph1 + dph3*jguess
228            mat4 = dph2 + dph3*iguess
229
230            determinant = mat1*mat4 - mat2*mat3
231
232            deli = (dthp*mat4 - mat2*dphp)/determinant
233            delj = (mat1*dphp - dthp*mat3)/determinant
234
235            if (abs(deli) < converge .and. 
236     &          abs(delj) < converge) exit iter_loop1
237
238            iguess = iguess + deli
239            jguess = jguess + delj
240
241          end do iter_loop1
242
243          if (iter <= max_iter) then
244
245            !***
246            !*** successfully found i,j - compute weights
247            !***
248
249            wgts(1) = (one-iguess)*(one-jguess)
250            wgts(2) = iguess*(one-jguess)
251            wgts(3) = iguess*jguess
252            wgts(4) = (one-iguess)*jguess
253            call store_link_bilin(dst_add, src_add, wgts, nmap)
254          else
255            write(nulou,*) 'Point coords: ',plat,plon
256            write(nulou,*) 'Dest grid lats: ',src_lats
257            write(nulou,*) 'Dest grid lons: ',src_lons
258            write(nulou,*) 'Dest grid addresses: ',src_add
259            write(nulou,*) 'Current i,j : ',iguess, jguess
260            write(nulou,*) 
261     &          'Iteration for i,j exceed max iteration count'
262            stop
263          endif
264
265        else if (src_add(1) < 0) THEN
266
267          !***
268          !*** We are in the first or last bin or at least one of the 4
269          !*** neighbours was masked. Do distance-weighted average using
270          !*** the non-masked points among the 4 closest ones.
271
272          IF (nlogprt .eq. 2) then
273              WRITE(nulou,*) ' '
274              WRITE(nulou,*) 
275     &   'WARNING: Cannot make bilinear interpolation for target point'
276     &            ,dst_add
277              WRITE(nulou,*) 
278     &    'Using non-masked points among 4 nearest neighbors.'
279              CALL OASIS_FLUSH_SCRIP(nulou)
280          ENDIF
281
282          !***
283          !*** Find the 4 closest points
284          !***
285
286          coslat_dst = cos(plat)
287          sinlat_dst = sin(plat)
288          coslon_dst = cos(plon)
289          sinlon_dst = sin(plon)
290          src_add = 0
291          dist_min = bignum
292          src_lats = bignum
293          IF (min_add == 0) min_add = 1
294          IF (max_add > grid1_size) max_add = grid1_size
295          do srch_add = min_add,max_add
296            arg = coslat_dst*cos(grid1_center_lat(srch_add))*
297     &           (coslon_dst*cos(grid1_center_lon(srch_add)) +
298     &            sinlon_dst*sin(grid1_center_lon(srch_add)))+
299     &            sinlat_dst*sin(grid1_center_lat(srch_add))
300            IF (arg < -1.0d0) THEN
301                arg = -1.0d0
302            ELSE IF (arg > 1.0d0) THEN
303                arg = 1.0d0
304            END IF   
305            distance=acos(arg)
306
307            if (distance < dist_min) then
308                sort_loop: do n=1,4
309                if (distance < src_lats(n)) then
310                    do i=4,n+1,-1
311                      src_add (i) = src_add (i-1)
312                      src_lats(i) = src_lats(i-1)
313                    end do
314                    src_add (n) = srch_add
315                    src_lats(n) = distance
316                    dist_min = src_lats(4)
317                    exit sort_loop
318                endif
319                end do sort_loop
320            endif
321          end do
322
323          src_lons = one/(src_lats + tiny)
324          distance = sum(src_lons)
325          src_lats = src_lons/distance
326
327          !***
328          !*** Among 4 closest points, keep only the non-masked ones
329          !***
330
331          icount = 0
332          do n=1,4
333            if (grid1_mask(src_add(n)) .or.
334     &          (.not. grid1_mask(src_add(n)) .and. lextrapdone)) then
335                icount = icount + 1
336            else
337                src_lats(n) = zero
338            endif
339          end do
340
341          if (icount > 0) then
342              !*** renormalize weights
343              sum_wgts = sum(src_lats)
344              wgts(1) = src_lats(1)/sum_wgts
345              wgts(2) = src_lats(2)/sum_wgts
346              wgts(3) = src_lats(3)/sum_wgts
347              wgts(4) = src_lats(4)/sum_wgts
348
349              grid2_frac(dst_add) = one
350              call store_link_bilin(dst_add, src_add, wgts, nmap)
351          ELSE
352              IF (ll_nnei .eqv. .true. ) then
353              IF (nlogprt .ge. 2) THEN
354                  WRITE(nulou,*) '  '
355                  WRITE(nulou,*) 
356     &                'All 4 surrounding source grid points are masked'
357                  WRITE(nulou,*) 'for target point ',dst_add
358                  WRITE(nulou,*) 'with longitude and latitude', 
359     &                plon, plat
360                  WRITE(nulou,*) 
361     &                'Using the nearest non-masked neighbour.'
362                  CALL OASIS_FLUSH_SCRIP(nulou)
363              ENDIF
364              src_latsnn = bignum
365!cdir novector
366              do srch_add = min_add,max_add
367                if (grid1_mask(srch_add) .or.
368     &          (.not. grid1_mask(srch_add) .and. lextrapdone)) THEN
369                    arg = coslat_dst*cos(grid1_center_lat(srch_add))*
370     &                  (coslon_dst*cos(grid1_center_lon(srch_add)) +
371     &                  sinlon_dst*sin(grid1_center_lon(srch_add)))+
372     &                  sinlat_dst*sin(grid1_center_lat(srch_add))
373                    IF (arg < -1.0d0) THEN
374                        arg = -1.0d0
375                    ELSE IF (arg > 1.0d0) THEN
376                        arg = 1.0d0
377                    END IF
378                    distance=acos(arg)
379
380                    if (distance < src_latsnn) then
381                        src_addnn = srch_add
382                        src_latsnn = distance
383                    endif
384                endif
385              end DO
386              IF (nlogprt .ge. 2) THEN
387                  WRITE(nulou,*) 
388     &                'Nearest non masked neighbour is source point '
389     &                ,src_addnn
390                  WRITE(nulou,*) 'with longitude and latitude',
391     &                grid1_center_lon(src_addnn), 
392     &                grid1_center_lat(src_addnn) 
393                  WRITE(nulou,*) '  '
394                  CALL OASIS_FLUSH_SCRIP(nulou)
395              ENDIF
396              wgts(1) = 1.
397              wgts(2) = 0.
398              wgts(3) = 0.
399              wgts(4) = 0.
400              src_add(1) = src_addnn
401              src_add(2) = 0
402              src_add(3) = 0
403              src_add(4) = 0
404
405              grid2_frac(dst_add) = one
406              call store_link_bilin(dst_add, src_add, wgts, nmap)
407          endif
408          ENDIF
409      ENDIF
410      end do grid_loop1
411!
412!-----------------------------------------------------------------------
413
414      end subroutine remap_bilin_reduced
415
416!***********************************************************************
417
418      subroutine grid_search_bilin_rd(src_add, src_lats, src_lons,
419     &                                min_add, max_add,
420     &                                plat, plon, src_grid_dims,
421     &                                src_center_lat, src_center_lon,
422     &                                lextrapdone)
423
424!-----------------------------------------------------------------------
425!
426!     this routine finds the location of the search point plat, plon
427!     in the source grid and returns the corners needed for a bilinear
428!     interpolation. The target grid is a reduced grid.
429!
430!-----------------------------------------------------------------------
431
432!-----------------------------------------------------------------------
433!
434!     output variables
435!
436!-----------------------------------------------------------------------
437
438      integer (kind=int_kind), dimension(4), intent(out) ::
439     &        src_add  ! address of each corner point enclosing P
440
441      real (kind=dbl_kind), dimension(4), intent(out) ::
442     &        src_lats, ! latitudes  of the four corner points
443     &        src_lons  ! longitudes of the four corner points
444
445      integer (kind=int_kind) :: min_add, max_add
446
447!-----------------------------------------------------------------------
448!
449!     input variables
450!
451!-----------------------------------------------------------------------
452
453      real (kind=dbl_kind), intent(in) ::
454     &        plat,   ! latitude  of the search point
455     &        plon    ! longitude of the search point
456
457      integer (kind=int_kind), dimension(2), intent(in) ::
458     &        src_grid_dims  ! size of each src grid dimension
459
460      real (kind=dbl_kind), dimension(:), intent(in) ::
461     &        src_center_lat, ! latitude  of each src grid center
462     &        src_center_lon  ! longitude of each src grid center
463
464      LOGICAL ::
465     &           lextrapdone   ! logical, true if EXTRAP done on field
466
467!-----------------------------------------------------------------------
468!
469!     local variables
470!
471!-----------------------------------------------------------------------
472
473      integer (kind=int_kind) :: n, next_n, srch_add, ni,  ! dummy indices
474     &    nx, ny, ntotmask,           ! dimensions of src grid
475     &    inter_add ! add for restricting search
476!
477      integer (kind=int_kind), DIMENSION(4) :: src_bid
478
479!-----------------------------------------------------------------------
480!
481!     restrict search first using bins
482!
483!-----------------------------------------------------------------------
484
485      nx = src_grid_dims(1)
486      inter_add = 0
487
488      src_add = 0
489
490      min_add = size(src_center_lat) + 1
491      max_add = 1
492      if (plat >= bin_lats_r(1,1)) then
493          min_add = 0
494          max_add = bin_addr1_r(4,1)
495          inter_add = bin_addr1_r(3,1)
496      else if (plat <= bin_lats_r(1,num_srch_red)) then
497          max_add = nx + 1
498          min_add = bin_addr1_r(1,num_srch_red)
499          inter_add = bin_addr1_r(3,num_srch_red)
500      else
501          srch_loopred: do n=1,num_srch_red
502            if (plat <= bin_lats_r(1,n) 
503     &        .and. plat >= bin_lats_r(2,n)) then
504                min_add = bin_addr1_r(1,n)
505                max_add = bin_addr1_r(4,n)
506                inter_add = bin_addr1_r(3,n)
507                exit srch_loopred
508            endif
509          end DO srch_loopred
510      ENDIF
511
512!-----------------------------------------------------------------------
513!
514!     now perform a more detailed search
515!
516!-----------------------------------------------------------------------
517      if (min_add .ne. 0 .and. max_add .ne. nx+1) THEN
518         !*** The concerned bins are not the top north or south ones.
519         !*** We should be able to find the four corners
520         !*** for the bilinear interpolation.
521
522          IF ( plon <= src_center_lon(min_add) ) THEN
523              src_add(1) = inter_add-1
524              src_add(2) = min_add
525          ELSE IF ( plon > src_center_lon(inter_add-1) ) then
526              src_add(1) = inter_add-1
527              src_add(2) = min_add
528          else
529              srch_loop2: do srch_add = min_add, inter_add-2
530                 if ( (plon > src_center_lon(srch_add)) .and.
531     &            (plon <= src_center_lon(srch_add+1)) )then
532                     src_add(1) = srch_add
533                     src_add(2) = srch_add+1
534                     exit srch_loop2
535                 endif
536               end do srch_loop2
537           ENDIF
538           IF ( plon <= src_center_lon(inter_add) ) THEN
539               src_add(3) = max_add
540               src_add(4) = inter_add
541           ELSE IF ( plon >= src_center_lon(max_add) ) then
542               src_add(3) = max_add
543               src_add(4) = inter_add
544           else
545               srch_loop3: do srch_add = inter_add, max_add
546                 if ( (plon >= src_center_lon(srch_add)) .and.
547     &             (plon <= src_center_lon(srch_add+1)) )then
548                      src_add(3) = srch_add
549                      src_add(4) = srch_add+1
550                      exit srch_loop3
551                  endif
552               enddo srch_loop3
553           ENDIF
554           src_lats(1) = src_center_lat(src_add(3))
555           src_lats(2) = src_center_lat(src_add(4))
556           src_lats(3) = src_center_lat(src_add(2))
557           src_lats(4) = src_center_lat(src_add(1))
558     
559           src_lons(1) = src_center_lon(src_add(3))
560           src_lons(2) = src_center_lon(src_add(4))
561           src_lons(3) = src_center_lon(src_add(2))
562           src_lons(4) = src_center_lon(src_add(1))
563
564           src_bid=src_add
565
566           src_add(1) = src_bid(3)
567           src_add(2) = src_bid(4)
568           src_add(3) = src_bid(2)
569           src_add(4) = src_bid(1)
570   
571           ! Check if one point is masked; IF so, nearest-neighbour
572           ! interpolation will be used
573
574           ntotmask = 0
575           do ni=1,4
576             if (.not. grid1_mask(src_add(ni)).and. 
577     &           .not. lextrapdone) 
578     &           ntotmask = ntotmask + 1 
579           end DO
580           IF (ntotmask .gt. 0) src_add(1) = -src_add(1) 
581         
582       ELSE 
583
584           !*** We are in the first or last bin.  Put src_add = -1 so that
585           !*** distance-weighted average of the 4 non-masked closest points
586           !*** is done in calling routine.
587
588           src_add = -1
589
590       ENDIF
591
592!-----------------------------------------------------------------------
593
594      end subroutine grid_search_bilin_rd 
595
596!***********************************************************************
597
598      subroutine store_link_bilin(dst_add, src_add, weights, nmap)
599
600!-----------------------------------------------------------------------
601!
602!     this routine stores the address and weight for four links
603!     associated with one destination point in the appropriate address
604!     and weight arrays and resizes those arrays if necessary.
605!
606!-----------------------------------------------------------------------
607
608!-----------------------------------------------------------------------
609!
610!     input variables
611!
612!-----------------------------------------------------------------------
613
614      integer (kind=int_kind), intent(in) ::
615     &        dst_add,  ! address on destination grid
616     &        nmap      ! identifies which direction for mapping
617
618      integer (kind=int_kind), dimension(4), intent(in) ::
619     &        src_add   ! addresses on source grid
620
621      real (kind=dbl_kind), dimension(4), intent(in) ::
622     &        weights ! array of remapping weights for these links
623
624!-----------------------------------------------------------------------
625!
626!     local variables
627!
628!-----------------------------------------------------------------------
629
630      integer (kind=int_kind) :: n, ! dummy index
631     &       num_links_old          ! placeholder for old link number
632
633!-----------------------------------------------------------------------
634!
635!     increment number of links and check to see if remap arrays need
636!     to be increased to accomodate the new link.  then store the
637!     link.
638!
639!-----------------------------------------------------------------------
640
641      select case (nmap)
642      case(1)
643
644        num_links_old  = num_links_map1
645        num_links_map1 = num_links_old + 4
646
647        if (num_links_map1 > max_links_map1) 
648     &     call resize_remap_vars(1,resize_increment)
649
650        do n=1,4
651          grid1_add_map1(num_links_old+n) = src_add(n)
652          grid2_add_map1(num_links_old+n) = dst_add
653          wts_map1    (1,num_links_old+n) = weights(n)
654        end do
655
656      case(2)
657
658        num_links_old  = num_links_map2
659        num_links_map2 = num_links_old + 4
660
661        if (num_links_map2 > max_links_map2) 
662     &     call resize_remap_vars(2,resize_increment)
663
664        do n=1,4
665          grid1_add_map2(num_links_old+n) = dst_add
666          grid2_add_map2(num_links_old+n) = src_add(n)
667          wts_map2    (1,num_links_old+n) = weights(n)
668        end do
669
670      end select
671
672!-----------------------------------------------------------------------
673
674      end subroutine store_link_bilin
675
676!***********************************************************************
677
678      end module remap_bilinear_reduced
679
680!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Note: See TracBrowser for help on using the repository browser.