source: CPL/oasis3-mct/branches/OASIS3-MCT_2.0_branch/lib/scrip/src/scrip.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: 18.0 KB
Line 
1C****
2C               *****************************
3C               * OASIS ROUTINE  -  LEVEL 1 *
4C               * -------------     ------- *
5C               *****************************
6C****
7C***********************************************************************
8C     This routine belongs to the SCRIP library. It is modified to run
9C     within OASIS.
10C     Modifications:
11C       - routine does not read namelist but gets parameters from the
12C         calling routine scriprmp.f
13C       - map-method and noralize-option are written in capital letters
14C       - routine grid_init is not called from scrip any more but was
15C         called earlier from scriprmp
16C       - call of two extra routines: free_grids and free_remap_vars to 
17C         allow multiple calls of SCRIP
18C       - added case for GAUSWGT 
19C       - added 'REDUCED' case for bilinear and bicubic.
20C       - hard coded num_maps=1 for USE in OASIS
21C       - added lextrapdone argument
22C
23C     Modified by            V. Gayler,  M&D                  20.09.2001
24C     Modified by            D. Declat,  CERFACS              27.06.2002
25C     Modified by            S. Valcke,  CERFACS              27.08.2002
26C***********************************************************************
27!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
28!
29!     This routine is the driver for computing the addresses and weights
30!     for interpolating between two grids on a sphere.
31!
32!-----------------------------------------------------------------------
33!
34!     CVS:$Id: scrip.f 1831 2009-01-09 17:19:08Z valcke $
35!
36!     Copyright (c) 1997, 1998 the Regents of the University of
37!       California.
38!
39!     This software and ancillary information (herein called software)
40!     called SCRIP is made available under the terms described here. 
41!     The software has been approved for release with associated
42!     LA-CC Number 98-45.
43!
44!     Unless otherwise indicated, this software has been authored
45!     by an employee or employees of the University of California,
46!     operator of the Los Alamos National Laboratory under Contract
47!     No. W-7405-ENG-36 with the U.S. Department of Energy.  The U.S.
48!     Government has rights to use, reproduce, and distribute this
49!     software.  The public may copy and use this software without
50!     charge, provided that this Notice and any statement of authorship
51!     are reproduced on all copies.  Neither the Government nor the
52!     University makes any warranty, express or implied, or assumes
53!     any liability or responsibility for the use of this software.
54!
55!     If software is modified to produce derivative works, such modified
56!     software should be clearly marked, so as not to confuse it with
57!     the version available from Los Alamos National Laboratory.
58!
59!***********************************************************************
60
61      subroutine scrip
62     $           (interp_file1, map1_name, m_method, n_opt, 
63     $            lextrapdone, rl_varmul, id_scripvoi)
64
65!-----------------------------------------------------------------------
66
67      use kinds_mod                  ! module defining data types
68      use constants                  ! module for common constants
69      use iounits                    ! I/O unit manager
70      use timers                     ! CPU timers
71      use grids                      ! module with grid information
72      use remap_vars                 ! common remapping variables
73      use remap_conservative         ! routines for conservative remap
74      use remap_distance_weight      ! routines for dist-weight remap
75      use remap_gaussian_weight      ! routines for gaus-weight remap
76      use remap_bilinear             ! routines for bilinear interp
77      use remap_bicubic              ! routines for bicubic  interp
78      use remap_bilinear_reduced     ! routines for bilinear interp
79      use remap_bicubic_reduced      ! routines for bicubic interp
80      use remap_write                ! routines for remap output
81      USE mod_oasis_flush
82
83      implicit none
84
85!-----------------------------------------------------------------------
86!
87!     input variables
88!
89!-----------------------------------------------------------------------
90
91      character (char_len), intent(in) ::
92     &           interp_file1, ! filename for output remap data (map1)
93     &           map1_name     ! name for mapping from grid1 to grid2
94
95      character*8, intent(in) ::
96     &           m_method,     ! choice for mapping method
97     &           n_opt         ! option for normalizing weights
98
99      LOGICAL ::
100     &           lextrapdone   ! logical, true if EXTRAP done on field
101
102      REAL (kind=real_kind) ::
103     &    rl_varmul             ! Gaussian variance (for GAUSWGT)
104
105      INTEGER (kind=int_kind) ::
106     &     id_scripvoi          ! number of neighbours for DISTWGT and GAUSWGT
107
108!-----------------------------------------------------------------------
109!
110!     local variables
111!
112!-----------------------------------------------------------------------
113
114      integer (kind=int_kind) ::
115     &           n             ! dummy counter
116
117      character (char_len) ::
118     &           interp_file2, ! filename for output remap data (map2)
119     &           map2_name,    ! name for mapping from grid2 to grid1
120     &           output_opt,   ! option for output conventions
121     &           map_method,   ! choice for mapping method
122     &           normalize_opt ! option for normalizing weights
123!-----------------------------------------------------------------------
124!
125      IF (nlogprt .GE. 2) THEN
126         WRITE (UNIT = nulou,FMT = *)' '
127         WRITE (UNIT = nulou,FMT = *)'Entering routine scrip'
128         CALL OASIS_FLUSH_SCRIP(nulou)
129      ENDIF
130!
131!-----------------------------------------------------------------------
132!
133!     initialize timers
134!
135!-----------------------------------------------------------------------
136
137      call timers_init
138      do n=1,max_timers
139        call timer_clear(n)
140      end do
141
142!-----------------------------------------------------------------------
143!
144!     initialize variables of former SCRIP namelist
145!
146!-----------------------------------------------------------------------
147
148      interp_file2  = 'unknown'
149      map2_name     = 'unknown'
150      luse_grid1_area = .false.
151      luse_grid2_area = .false.
152      num_maps      = 1
153      output_opt    = 'scrip'
154
155      map_method = m_method
156      normalize_opt = n_opt
157
158      select case(map_method)
159      case ('CONSERV')
160        map_type = map_type_conserv
161      case ('BILINEAR')
162        map_type = map_type_bilinear
163      case ('BICUBIC')
164        map_type = map_type_bicubic
165      case ('DISTWGT')
166        map_type = map_type_distwgt
167      case ('GAUSWGT')
168        map_type = map_type_gauswgt
169      case default
170        stop 'unknown mapping method'
171      end select
172
173      SELECT CASE (normalize_opt)
174      CASE ('FRACNNEI')
175          lfracnnei = .true.
176      END SELECT
177
178      select case(normalize_opt(1:4))
179      case ('NONE')
180        norm_opt = norm_opt_none
181      case ('FRAC')
182        norm_opt = norm_opt_frcarea
183      case ('DEST')
184        norm_opt = norm_opt_dstarea
185      CASE ('NONO')
186        norm_opt = norm_opt_nonorm
187      case default
188        stop 'unknown normalization option'
189      end select
190!
191      IF (nlogprt .GE. 2) THEN
192         WRITE (UNIT = nulou,FMT = *)' Computing remappings between: '
193     &        ,grid1_name
194         WRITE (UNIT = nulou,FMT = *)'                          and  '
195     &       ,grid2_name
196         CALL OASIS_FLUSH_SCRIP(nulou)
197      ENDIF
198!
199!-----------------------------------------------------------------------
200!
201!     initialize some remapping variables.
202!
203!-----------------------------------------------------------------------
204      call init_remap_vars (id_scripvoi)
205!-----------------------------------------------------------------------
206!
207!     call appropriate interpolation setup routine based on type of
208!     remapping requested.
209!
210!-----------------------------------------------------------------------
211      select case(map_type)
212      case(map_type_conserv)
213          call remap_conserv
214      case(map_type_bilinear)
215          IF (restrict_TYPE == 'REDUCED') then
216              call remap_bilin_reduced (lextrapdone)
217          ELSE
218              call remap_bilin (lextrapdone)
219          endif
220      case(map_type_distwgt)
221          call remap_distwgt (lextrapdone, id_scripvoi)
222      case(map_type_gauswgt)
223          call remap_gauswgt (lextrapdone, id_scripvoi, rl_varmul)
224      case(map_type_bicubic)
225          IF (restrict_TYPE == 'REDUCED') then
226              call remap_bicub_reduced (lextrapdone)
227          ELSE
228              call remap_bicub (lextrapdone)
229          endif
230       case default
231           stop 'Invalid Map Type'
232      end select
233     
234      CALL sort_add(grid2_add_map1, grid1_add_map1, wts_map1,
235     $   num_links_map1, num_wts)
236      IF (lfracnnei) THEN
237          CALL fracnnei(grid1_size, grid2_size, grid1_mask, grid2_mask,
238     $        grid1_center_lon, grid1_center_lat,
239     $        grid2_center_lon, grid2_center_lat,
240     $        num_links_map1, num_wts,
241     $        wts_map1, grid1_add_map1, grid2_add_map1)
242          lfracnnei = .false.
243      ENDIF
244!
245#ifdef TREAT_OVERLAY
246!
247! Change address if overlap point were found
248      IF (map_type == 1) THEN
249          DO n = 1, num_links_map1
250            IF (grid1_add_map1(n) .ne. 0) THEN
251                grid1_add_map1(n) = grid1_add_repl1(grid1_add_map1(n))
252            ENDIF
253          END DO
254      ENDIF
255!
256#endif
257!
258!
259!-----------------------------------------------------------------------
260!
261!     reduce size of remapping arrays and then write remapping info
262!     to a file.
263!
264!-----------------------------------------------------------------------
265      if (num_links_map1 /= max_links_map1) then
266        call resize_remap_vars(1, num_links_map1-max_links_map1)
267      endif
268      if ((num_maps > 1) .and. (num_links_map2 /= max_links_map2)) then
269        call resize_remap_vars(2, num_links_map2-max_links_map2)
270      endif
271
272      call write_remap(map1_name, map2_name, 
273     &                 interp_file1, interp_file2, output_opt)
274
275!-----------------------------------------------------------------------
276!
277!     deallocate allocatable arrays
278!
279!-----------------------------------------------------------------------
280
281      call free_grids
282      call free_remap_vars
283!
284      IF (nlogprt .GE. 2) THEN
285         WRITE (UNIT = nulou,FMT = *)' '
286         WRITE (UNIT = nulou,FMT = *)'Leaving routine scrip'
287         CALL OASIS_FLUSH_SCRIP(nulou)
288      ENDIF
289!-----------------------------------------------------------------------!
290      end subroutine scrip
291!
292      subroutine sort_add(add1, add2, weights, num_links, num_wts)
293
294!-----------------------------------------------------------------------
295!
296!     this routine sorts address and weight arrays based on the
297!     destination address with the source address as a secondary
298!     sorting criterion.  the method is a standard heap sort.
299!
300!-----------------------------------------------------------------------
301
302      use kinds_mod     ! defines common data types
303      use constants     ! defines common scalar constants
304      USE mod_oasis_flush
305
306      implicit none
307
308!-----------------------------------------------------------------------
309!
310!     Input and Output arrays
311!
312!-----------------------------------------------------------------------
313
314      integer (kind=int_kind), intent(in) :: num_links, num_wts
315      integer (kind=int_kind), intent(inout), dimension(num_links) ::
316     &        add1,       ! destination address array (num_links)
317     &        add2        ! source      address array
318
319      real (kind=dbl_kind), intent(inout),
320     &    dimension(num_wts, num_links) ::
321     &        weights     ! remapping weights (num_wts, num_links)
322
323
324!-----------------------------------------------------------------------
325!
326!     local variables
327!
328!-----------------------------------------------------------------------
329
330      integer (kind=int_kind) ::
331!     &          num_links,          ! num of links for this mapping
332!     &          num_wts,            ! num of weights for this mapping
333     &          add1_tmp, add2_tmp, ! temp for addresses during swap
334     &          nwgt,
335     &          lvl, final_lvl,     ! level indexes for heap sort levels
336     &          chk_lvl1, chk_lvl2, max_lvl
337
338      real (kind=dbl_kind), dimension(SIZE(weights,DIM=1)) ::
339     &          wgttmp              ! temp for holding wts during swap
340!-----------------------------------------------------------------------
341!
342      IF (nlogprt .GE. 2) THEN
343         WRITE (UNIT = nulou,FMT = *)' '
344         WRITE (UNIT = nulou,FMT = *)'Entering routine sort_add'
345         CALL OASIS_FLUSH_SCRIP(nulou)
346      ENDIF
347!
348!-----------------------------------------------------------------------
349!
350!     determine total number of links to sort and number of weights
351!
352!-----------------------------------------------------------------------
353
354!      num_links = SIZE(add1)
355!      num_wts   = SIZE(weights, DIM=1)
356
357!-----------------------------------------------------------------------
358!
359!     start at the lowest level (N/2) of the tree and sift lower
360!     values to the bottom of the tree, promoting the larger numbers
361!
362!-----------------------------------------------------------------------
363
364      do lvl=num_links/2,1,-1
365
366        final_lvl = lvl
367        add1_tmp = add1(lvl)
368        add2_tmp = add2(lvl)
369        wgttmp(:) = weights(:,lvl)
370
371        !***
372        !*** loop until proper level is found for this link, or reach
373        !*** bottom
374        !***
375
376        sift_loop1: do
377
378          !***
379          !*** find the largest of the two daughters
380          !***
381
382          chk_lvl1 = 2*final_lvl
383          chk_lvl2 = 2*final_lvl+1
384          if (chk_lvl1 .EQ. num_links) chk_lvl2 = chk_lvl1
385
386          if ((add1(chk_lvl1) >  add1(chk_lvl2)) .OR.
387     &       ((add1(chk_lvl1) == add1(chk_lvl2)) .AND.
388     &        (add2(chk_lvl1) >  add2(chk_lvl2)))) then
389            max_lvl = chk_lvl1
390          else
391            max_lvl = chk_lvl2
392          endif
393
394          !***
395          !*** if the parent is greater than both daughters,
396          !*** the correct level has been found
397          !***
398
399          if ((add1_tmp .GT. add1(max_lvl)) .OR.
400     &       ((add1_tmp .EQ. add1(max_lvl)) .AND.
401     &        (add2_tmp .GT. add2(max_lvl)))) then
402            add1(final_lvl) = add1_tmp
403            add2(final_lvl) = add2_tmp
404            weights(:,final_lvl) = wgttmp(:)
405            exit sift_loop1
406
407          !***
408          !*** otherwise, promote the largest daughter and push
409          !*** down one level in the tree.  if haven't reached
410          !*** the end of the tree, repeat the process.  otherwise
411          !*** store last values and exit the loop
412          !***
413
414          else
415            add1(final_lvl) = add1(max_lvl)
416            add2(final_lvl) = add2(max_lvl)
417            weights(:,final_lvl) = weights(:,max_lvl)
418
419            final_lvl = max_lvl
420            if (2*final_lvl > num_links) then
421              add1(final_lvl) = add1_tmp
422              add2(final_lvl) = add2_tmp
423              weights(:,final_lvl) = wgttmp(:)
424              exit sift_loop1
425            endif
426          endif
427        end do sift_loop1
428      end do
429
430!-----------------------------------------------------------------------
431!
432!     now that the heap has been sorted, strip off the top (largest)
433!     value and promote the values below
434!
435!-----------------------------------------------------------------------
436
437      do lvl=num_links,3,-1
438
439        !***
440        !*** move the top value and insert it into the correct place
441        !***
442
443        add1_tmp = add1(lvl)
444        add1(lvl) = add1(1)
445
446        add2_tmp = add2(lvl)
447        add2(lvl) = add2(1)
448
449        wgttmp(:) = weights(:,lvl)
450        weights(:,lvl) = weights(:,1)
451
452        !***
453        !*** as above this loop sifts the tmp values down until proper
454        !*** level is reached
455        !***
456
457        final_lvl = 1
458
459        sift_loop2: do
460
461          !***
462          !*** find the largest of the two daughters
463          !***
464
465          chk_lvl1 = 2*final_lvl
466          chk_lvl2 = 2*final_lvl+1
467          if (chk_lvl2 >= lvl) chk_lvl2 = chk_lvl1
468
469          if ((add1(chk_lvl1) >  add1(chk_lvl2)) .OR.
470     &       ((add1(chk_lvl1) == add1(chk_lvl2)) .AND.
471     &        (add2(chk_lvl1) >  add2(chk_lvl2)))) then
472            max_lvl = chk_lvl1
473          else
474            max_lvl = chk_lvl2
475          endif
476
477          !***
478          !*** if the parent is greater than both daughters,
479          !*** the correct level has been found
480          !***
481
482          if ((add1_tmp >  add1(max_lvl)) .OR.
483     &       ((add1_tmp == add1(max_lvl)) .AND.
484     &        (add2_tmp >  add2(max_lvl)))) then
485            add1(final_lvl) = add1_tmp
486            add2(final_lvl) = add2_tmp
487            weights(:,final_lvl) = wgttmp(:)
488            exit sift_loop2
489
490          !***
491          !*** otherwise, promote the largest daughter and push
492          !*** down one level in the tree.  if haven't reached
493          !*** the end of the tree, repeat the process.  otherwise
494          !*** store last values and exit the loop
495          !***
496
497          else
498            add1(final_lvl) = add1(max_lvl)
499            add2(final_lvl) = add2(max_lvl)
500            weights(:,final_lvl) = weights(:,max_lvl)
501
502            final_lvl = max_lvl
503            if (2*final_lvl >= lvl) then
504              add1(final_lvl) = add1_tmp
505              add2(final_lvl) = add2_tmp
506              weights(:,final_lvl) = wgttmp(:)
507              exit sift_loop2
508            endif
509          endif
510        end do sift_loop2
511      end do
512
513      !***
514      !*** swap the last two entries
515      !***
516
517
518      add1_tmp = add1(2)
519      add1(2)  = add1(1)
520      add1(1)  = add1_tmp
521
522      add2_tmp = add2(2)
523      add2(2)  = add2(1)
524      add2(1)  = add2_tmp
525
526      wgttmp (:)   = weights(:,2)
527      weights(:,2) = weights(:,1)
528      weights(:,1) = wgttmp (:)
529!
530      IF (nlogprt .GE. 2) THEN
531         WRITE (UNIT = nulou,FMT = *)' '
532         WRITE (UNIT = nulou,FMT = *)'Leaving routine sort_add'
533         CALL OASIS_FLUSH_SCRIP(nulou)
534      ENDIF
535!
536!-----------------------------------------------------------------------
537
538      end subroutine sort_add
539
540!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Note: See TracBrowser for help on using the repository browser.