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

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

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

File size: 16.5 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,v 1.1.1.1 2005/03/23 16:01:12 adm Exp $
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
82      implicit none
83
84!-----------------------------------------------------------------------
85!
86!     input variables
87!
88!-----------------------------------------------------------------------
89
90      character (char_len), intent(in) ::
91     &           interp_file1, ! filename for output remap data (map1)
92     &           map1_name     ! name for mapping from grid1 to grid2
93
94      character*8, intent(in) ::
95     &           m_method,     ! choice for mapping method
96     &           n_opt         ! option for normalizing weights
97
98      LOGICAL ::
99     &           lextrapdone   ! logical, true if EXTRAP done on field
100
101      REAL (kind=real_kind) ::
102     &    rl_varmul             ! Gaussian variance (for GAUSWGT)
103
104      INTEGER (kind=int_kind) ::
105     &     id_scripvoi          ! number of neighbours for DISTWGT and GAUSWGT
106
107!-----------------------------------------------------------------------
108!
109!     local variables
110!
111!-----------------------------------------------------------------------
112
113      integer (kind=int_kind) ::
114     &           n             ! dummy counter
115
116      character (char_len) ::
117     &           interp_file2, ! filename for output remap data (map2)
118     &           map2_name,    ! name for mapping from grid2 to grid1
119     &           output_opt,   ! option for output conventions
120     &           map_method,   ! choice for mapping method
121     &           normalize_opt ! option for normalizing weights
122
123!-----------------------------------------------------------------------
124!
125!     initialize timers
126!
127!-----------------------------------------------------------------------
128
129      call timers_init
130      do n=1,max_timers
131        call timer_clear(n)
132      end do
133
134!-----------------------------------------------------------------------
135!
136!     initialize variables of former SCRIP namelist
137!
138!-----------------------------------------------------------------------
139
140      interp_file2  = 'unknown'
141      map2_name     = 'unknown'
142      luse_grid1_area = .false.
143      luse_grid2_area = .false.
144      num_maps      = 1
145      output_opt    = 'scrip'
146
147      map_method = m_method
148      normalize_opt = n_opt
149
150      select case(map_method)
151      case ('CONSERV')
152        map_type = map_type_conserv
153      case ('BILINEAR')
154        map_type = map_type_bilinear
155      case ('BICUBIC')
156        map_type = map_type_bicubic
157      case ('DISTWGT')
158        map_type = map_type_distwgt
159      case ('GAUSWGT')
160        map_type = map_type_gauswgt
161      case default
162        stop 'unknown mapping method'
163      end select
164
165      select case(normalize_opt(1:4))
166      case ('NONE')
167        norm_opt = norm_opt_none
168      case ('FRAC')
169        norm_opt = norm_opt_frcarea
170      case ('DEST')
171        norm_opt = norm_opt_dstarea
172      CASE ('NONO')
173        norm_opt = norm_opt_nonorm
174      case default
175        stop 'unknown normalization option'
176      end select
177
178      write(stdout, *) ' Computing remappings between: ',grid1_name
179      write(stdout, *) '                          and  ',grid2_name
180
181!-----------------------------------------------------------------------
182!
183!     initialize some remapping variables.
184!
185!-----------------------------------------------------------------------
186
187      call init_remap_vars (id_scripvoi)
188
189!-----------------------------------------------------------------------
190!
191!     call appropriate interpolation setup routine based on type of
192!     remapping requested.
193!
194!-----------------------------------------------------------------------
195
196      select case(map_type)
197      case(map_type_conserv)
198          call remap_conserv
199      case(map_type_bilinear)
200          IF (restrict_TYPE == 'REDUCED') then
201              call remap_bilin_reduced (lextrapdone)
202          ELSE
203              call remap_bilin (lextrapdone)
204          endif
205      case(map_type_distwgt)
206          call remap_distwgt (lextrapdone, id_scripvoi)
207      case(map_type_gauswgt)
208          call remap_gauswgt (lextrapdone, id_scripvoi, rl_varmul)
209      case(map_type_bicubic)
210          IF (restrict_TYPE == 'REDUCED') then
211              call remap_bicub_reduced (lextrapdone)
212          ELSE
213              call remap_bicub (lextrapdone)
214          endif
215       case default
216           stop 'Invalid Map Type'
217      end select
218      CALL sort_add(grid2_add_map1, grid1_add_map1, wts_map1,
219     $   num_links_map1, num_wts)
220      IF (lfracnnei) THEN
221          CALL fracnnei(grid1_size, grid2_size, grid1_mask, grid2_mask,
222     $        grid1_center_lon, grid1_center_lat,
223     $        grid2_center_lon, grid2_center_lat,
224     $        num_links_map1, num_wts,
225     $        wts_map1, grid1_add_map1, grid2_add_map1)
226          lfracnnei = .false.
227      ENDIF
228
229!-----------------------------------------------------------------------
230!
231!     reduce size of remapping arrays and then write remapping info
232!     to a file.
233!
234!-----------------------------------------------------------------------
235      if (num_links_map1 /= max_links_map1) then
236        call resize_remap_vars(1, num_links_map1-max_links_map1)
237      endif
238      if ((num_maps > 1) .and. (num_links_map2 /= max_links_map2)) then
239        call resize_remap_vars(2, num_links_map2-max_links_map2)
240      endif
241      call write_remap(map1_name, map2_name, 
242     &                 interp_file1, interp_file2, output_opt)
243
244!-----------------------------------------------------------------------
245!
246!     deallocate allocatable arrays
247!
248!-----------------------------------------------------------------------
249
250      call free_grids
251      call free_remap_vars
252
253!-----------------------------------------------------------------------
254
255      end subroutine scrip
256!
257      subroutine sort_add(add1, add2, weights, num_links, num_wts)
258
259!-----------------------------------------------------------------------
260!
261!     this routine sorts address and weight arrays based on the
262!     destination address with the source address as a secondary
263!     sorting criterion.  the method is a standard heap sort.
264!
265!-----------------------------------------------------------------------
266
267      use kinds_mod     ! defines common data types
268      use constants     ! defines common scalar constants
269
270
271      implicit none
272
273!-----------------------------------------------------------------------
274!
275!     Input and Output arrays
276!
277!-----------------------------------------------------------------------
278
279      integer (kind=int_kind), intent(in) :: num_links, num_wts
280      integer (kind=int_kind), intent(inout), dimension(num_links) ::
281     &        add1,       ! destination address array (num_links)
282     &        add2        ! source      address array
283
284      real (kind=dbl_kind), intent(inout),
285     &    dimension(num_wts, num_links) ::
286     &        weights     ! remapping weights (num_wts, num_links)
287
288
289!-----------------------------------------------------------------------
290!
291!     local variables
292!
293!-----------------------------------------------------------------------
294
295      integer (kind=int_kind) ::
296!     &          num_links,          ! num of links for this mapping
297!     &          num_wts,            ! num of weights for this mapping
298     &          add1_tmp, add2_tmp, ! temp for addresses during swap
299     &          nwgt,
300     &          lvl, final_lvl,     ! level indexes for heap sort levels
301     &          chk_lvl1, chk_lvl2, max_lvl
302
303      real (kind=dbl_kind), dimension(SIZE(weights,DIM=1)) ::
304     &          wgttmp              ! temp for holding wts during swap
305
306!-----------------------------------------------------------------------
307!
308!     determine total number of links to sort and number of weights
309!
310!-----------------------------------------------------------------------
311
312!      num_links = SIZE(add1)
313!      num_wts   = SIZE(weights, DIM=1)
314
315!-----------------------------------------------------------------------
316!
317!     start at the lowest level (N/2) of the tree and sift lower
318!     values to the bottom of the tree, promoting the larger numbers
319!
320!-----------------------------------------------------------------------
321
322      do lvl=num_links/2,1,-1
323
324        final_lvl = lvl
325        add1_tmp = add1(lvl)
326        add2_tmp = add2(lvl)
327        wgttmp(:) = weights(:,lvl)
328
329        !***
330        !*** loop until proper level is found for this link, or reach
331        !*** bottom
332        !***
333
334        sift_loop1: do
335
336          !***
337          !*** find the largest of the two daughters
338          !***
339
340          chk_lvl1 = 2*final_lvl
341          chk_lvl2 = 2*final_lvl+1
342          if (chk_lvl1 .EQ. num_links) chk_lvl2 = chk_lvl1
343
344          if ((add1(chk_lvl1) >  add1(chk_lvl2)) .OR.
345     &       ((add1(chk_lvl1) == add1(chk_lvl2)) .AND.
346     &        (add2(chk_lvl1) >  add2(chk_lvl2)))) then
347            max_lvl = chk_lvl1
348          else
349            max_lvl = chk_lvl2
350          endif
351
352          !***
353          !*** if the parent is greater than both daughters,
354          !*** the correct level has been found
355          !***
356
357          if ((add1_tmp .GT. add1(max_lvl)) .OR.
358     &       ((add1_tmp .EQ. add1(max_lvl)) .AND.
359     &        (add2_tmp .GT. add2(max_lvl)))) then
360            add1(final_lvl) = add1_tmp
361            add2(final_lvl) = add2_tmp
362            weights(:,final_lvl) = wgttmp(:)
363            exit sift_loop1
364
365          !***
366          !*** otherwise, promote the largest daughter and push
367          !*** down one level in the tree.  if haven't reached
368          !*** the end of the tree, repeat the process.  otherwise
369          !*** store last values and exit the loop
370          !***
371
372          else
373            add1(final_lvl) = add1(max_lvl)
374            add2(final_lvl) = add2(max_lvl)
375            weights(:,final_lvl) = weights(:,max_lvl)
376
377            final_lvl = max_lvl
378            if (2*final_lvl > num_links) then
379              add1(final_lvl) = add1_tmp
380              add2(final_lvl) = add2_tmp
381              weights(:,final_lvl) = wgttmp(:)
382              exit sift_loop1
383            endif
384          endif
385        end do sift_loop1
386      end do
387
388!-----------------------------------------------------------------------
389!
390!     now that the heap has been sorted, strip off the top (largest)
391!     value and promote the values below
392!
393!-----------------------------------------------------------------------
394
395      do lvl=num_links,3,-1
396
397        !***
398        !*** move the top value and insert it into the correct place
399        !***
400
401        add1_tmp = add1(lvl)
402        add1(lvl) = add1(1)
403
404        add2_tmp = add2(lvl)
405        add2(lvl) = add2(1)
406
407        wgttmp(:) = weights(:,lvl)
408        weights(:,lvl) = weights(:,1)
409
410        !***
411        !*** as above this loop sifts the tmp values down until proper
412        !*** level is reached
413        !***
414
415        final_lvl = 1
416
417        sift_loop2: do
418
419          !***
420          !*** find the largest of the two daughters
421          !***
422
423          chk_lvl1 = 2*final_lvl
424          chk_lvl2 = 2*final_lvl+1
425          if (chk_lvl2 >= lvl) chk_lvl2 = chk_lvl1
426
427          if ((add1(chk_lvl1) >  add1(chk_lvl2)) .OR.
428     &       ((add1(chk_lvl1) == add1(chk_lvl2)) .AND.
429     &        (add2(chk_lvl1) >  add2(chk_lvl2)))) then
430            max_lvl = chk_lvl1
431          else
432            max_lvl = chk_lvl2
433          endif
434
435          !***
436          !*** if the parent is greater than both daughters,
437          !*** the correct level has been found
438          !***
439
440          if ((add1_tmp >  add1(max_lvl)) .OR.
441     &       ((add1_tmp == add1(max_lvl)) .AND.
442     &        (add2_tmp >  add2(max_lvl)))) then
443            add1(final_lvl) = add1_tmp
444            add2(final_lvl) = add2_tmp
445            weights(:,final_lvl) = wgttmp(:)
446            exit sift_loop2
447
448          !***
449          !*** otherwise, promote the largest daughter and push
450          !*** down one level in the tree.  if haven't reached
451          !*** the end of the tree, repeat the process.  otherwise
452          !*** store last values and exit the loop
453          !***
454
455          else
456            add1(final_lvl) = add1(max_lvl)
457            add2(final_lvl) = add2(max_lvl)
458            weights(:,final_lvl) = weights(:,max_lvl)
459
460            final_lvl = max_lvl
461            if (2*final_lvl >= lvl) then
462              add1(final_lvl) = add1_tmp
463              add2(final_lvl) = add2_tmp
464              weights(:,final_lvl) = wgttmp(:)
465              exit sift_loop2
466            endif
467          endif
468        end do sift_loop2
469      end do
470
471      !***
472      !*** swap the last two entries
473      !***
474
475
476      add1_tmp = add1(2)
477      add1(2)  = add1(1)
478      add1(1)  = add1_tmp
479
480      add2_tmp = add2(2)
481      add2(2)  = add2(1)
482      add2(1)  = add2_tmp
483
484      wgttmp (:)   = weights(:,2)
485      weights(:,2) = weights(:,1)
486      weights(:,1) = wgttmp (:)
487
488!-----------------------------------------------------------------------
489
490      end subroutine sort_add
491
492!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Note: See TracBrowser for help on using the repository browser.