source: CPL/oasis3-mct/branches/OASIS3-MCT_5.0_branch/lib/scrip/src/remap_vars.F90 @ 6331

Last change on this file since 6331 was 6331, checked in by aclsce, 17 months ago

Moved oasis-mct_5.0 in oasis3-mct/branches directory.
File size: 19.8 KB
Line 
1!****
2!                    ************************
3!                    *     OASIS MODULE     *
4!                    *     ------------     *
5!                    ************************
6!****
7!***********************************************************************
8!     This module belongs to the SCRIP library. It is modified to run
9!     within OASIS.
10!     Modifications:
11!       - introduction of logical flags to allow multiple calls of SCRIP
12!       - deallocation of arrays not needed any more
13!       - added CASE for GAUSWGT
14!
15!     Modified by            V. Gayler,  M&D                  20.09.2001
16!     Modified by            D. Declat,  CERFACS              27.06.2002
17!***********************************************************************
18!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
19!
20!     this module contains necessary variables for remapping between
21!     two grids.  also routines for resizing and initializing these
22!     variables.
23!
24!-----------------------------------------------------------------------
25!
26!     CVS:$Id: remap_vars.F 2826 2010-12-10 11:14:21Z valcke $
27!
28!     Copyright (c) 1997, 1998 the Regents of the University of
29!       California.
30!
31!     This software and ancillary information (herein called software)
32!     called SCRIP is made available under the terms described here. 
33!     The software has been approved for release with associated
34!     LA-CC Number 98-45.
35!
36!     Unless otherwise indicated, this software has been authored
37!     by an employee or employees of the University of California,
38!     operator of the Los Alamos National Laboratory under Contract
39!     No. W-7405-ENG-36 with the U.S. Department of Energy.  The U.S.
40!     Government has rights to use, reproduce, and distribute this
41!     software.  The public may copy and use this software without
42!     charge, provided that this Notice and any statement of authorship
43!     are reproduced on all copies.  Neither the Government nor the
44!     University makes any warranty, express or implied, or assumes
45!     any liability or responsibility for the use of this software.
46!
47!     If software is modified to produce derivative works, such modified
48!     software should be clearly marked, so as not to confuse it with
49!     the version available from Los Alamos National Laboratory.
50!
51!***********************************************************************
52
53      module remap_vars
54
55      use kinds_mod
56      use constants
57      use grids
58      USE mod_oasis_flush
59
60      implicit none
61
62!-----------------------------------------------------------------------
63!
64!     module variables
65!
66!-----------------------------------------------------------------------
67
68      integer (kind=int_kind), parameter :: norm_opt_none    = 1,  &
69                                            norm_opt_dstarea = 2,  &
70                                            norm_opt_frcarea = 3,  &
71                                            norm_opt_dstartr = 4,  &
72                                            norm_opt_frcartr = 5,  &
73                                            norm_opt_nonorm  = 6
74
75      integer (kind=int_kind), parameter :: map_type_conserv  = 1, &
76                                            map_type_bilinear = 2, &
77                                            map_type_bicubic  = 3, &
78                                            map_type_distwgt  = 4, &
79                                            map_type_gauswgt  = 5, &
80                                            map_type_loccwgt  = 6
81
82      integer (kind=int_kind), parameter :: norm_locc_uniform  = 1, &
83                                            norm_locc_distwgt  = 2, &
84                                            norm_locc_gauswgt  = 3
85
86      integer (kind=int_kind), save :: max_links_map1,  & ! current size of link arrays
87                                       num_links_map1,  & ! actual number of links for remapping
88                                       max_links_map2,  & ! current size of link arrays
89                                       num_links_map2,  & ! actual number of links for remapping
90                                       num_maps,        & ! num of remappings for this grid pair
91                                       num_wts,         & ! num of weights used in remapping
92                                       map_type,        & ! identifier for remapping method
93                                       norm_opt,        & ! option for normalization (conserv only)
94                                       conserve_opt,    & ! option for conservation first or second
95                                       resize_increment,& ! default amount to increase array size
96                                       norm_locc           ! locally conservative normalisation
97
98      integer (kind=int_kind), dimension(:), allocatable, save :: &
99           grid1_add_map1, & ! grid1 address for each link in mapping 1
100           grid2_add_map1, & ! grid2 address for each link in mapping 1
101           grid1_add_map2, & ! grid1 address for each link in mapping 2
102           grid2_add_map2    ! grid2 address for each link in mapping 2
103#ifdef TREAT_OVERLAY
104      integer (kind=int_kind), dimension(:), allocatable, save :: &
105           grid1_add_repl1 ! grid1 address to use after overlap calculation
106#endif
107      real (kind=dbl_kind), dimension(:,:), allocatable, save :: &
108           wts_map1, & ! map weights for each link (num_wts,max_links)
109           wts_map2    ! map weights for each link (num_wts,max_links)
110
111      logical (kind=log_kind), save :: lfracnnei = .false.
112 
113      logical (kind=log_kind), save :: first_conserv = .true. ! flag to
114                                ! indicate, whether scrip is called from
115                                ! oasis for the first time
116      logical (kind=log_kind), save :: first_call = .true. ! flag used in
117                                ! remap_conserve (store_link_cnsrv)
118
119      TYPE thread_remap
120
121         INTEGER (kind=int_kind) :: max_links
122         INTEGER (kind=int_kind) :: num_links
123         INTEGER (kind=int_kind) :: start_pos
124         INTEGER (kind=int_kind) :: nb_resize = 0
125         INTEGER (kind=int_kind), dimension(:), POINTER :: grid1_add
126         INTEGER (kind=int_kind), dimension(:), POINTER :: grid2_add
127         real (kind=dbl_kind), dimension(:,:), POINTER :: wts
128
129         CONTAINS
130
131         procedure :: resize => thread_resize_remap_vars
132
133      END TYPE
134
135      TYPE(thread_remap), dimension(:), allocatable, save :: sga_remap
136
137! mpi library include file
138#include <mpif.h>         
139
140!***********************************************************************
141
142      contains
143
144!***********************************************************************
145
146      subroutine init_remap_vars (id_scripvoi)
147
148!-----------------------------------------------------------------------
149!
150!     this routine initializes some variables and provides an initial
151!     allocation of arrays (fairly large so frequent resizing
152!     unnecessary).
153!
154!-----------------------------------------------------------------------
155!
156!     input variables
157!
158!-----------------------------------------------------------------------
159
160      INTEGER (kind=int_kind):: id_scripvoi ! number of neighbours for DISTWGT and GAUSWGT
161!-----------------------------------------------------------------------
162!
163      IF (nlogprt .GE. 2) THEN
164         WRITE (UNIT = nulou,FMT = *)' '
165         WRITE (UNIT = nulou,FMT = *)'Entering routine init_remap_vars'
166         CALL OASIS_FLUSH_SCRIP(nulou)
167      ENDIF
168!
169!-----------------------------------------------------------------------
170!
171!     determine the number of weights
172!
173!-----------------------------------------------------------------------
174
175      select case (map_type)
176      case(map_type_conserv)
177        num_wts = 3
178      case(map_type_bilinear)
179        num_wts = 1
180      case(map_type_bicubic)
181          IF (restrict_type == 'REDUCED') THEN
182              num_wts = 1
183          ELSE
184              num_wts = 4
185          ENDIF
186      case(map_type_distwgt)
187        num_wts = 1
188      case(map_type_gauswgt)
189        num_wts = 1
190      case(map_type_loccwgt)
191        num_wts = 1
192      end select
193
194!-----------------------------------------------------------------------
195!
196!     initialize num_links and set max_links to four times the largest
197!     of the destination grid sizes initially (can be changed later).
198!     set a default resize increment to increase the size of link
199!     arrays if the number of links exceeds the initial size
200!   
201!-----------------------------------------------------------------------
202
203      num_links_map1 = 0
204      select case (map_type)
205      case(map_type_conserv)
206          !EM could this be optimized ?
207          max_links_map1 = 4*grid2_size
208      case(map_type_bilinear)
209          max_links_map1 = 4*grid2_size
210      case(map_type_bicubic)
211          IF (restrict_type == 'REDUCED') THEN
212              max_links_map1 = 16*grid2_size
213          ELSE
214              max_links_map1 = 4*grid2_size
215          ENDIF
216      case(map_type_distwgt)
217          max_links_map1 = id_scripvoi*grid2_size
218      case(map_type_gauswgt)
219          max_links_map1 = id_scripvoi*grid2_size
220      case(map_type_loccwgt)
221          max_links_map1 = id_scripvoi*grid1_size
222      END select
223
224      if (num_maps > 1) then
225        num_links_map2 = 0
226        max_links_map1 = max(4*grid1_size,4*grid2_size)
227        max_links_map2 = max_links_map1
228      endif
229
230      resize_increment = 0.1*max(grid1_size,grid2_size)
231
232!-----------------------------------------------------------------------
233!
234!     allocate address and weight arrays for mapping 1
235!   
236!-----------------------------------------------------------------------
237
238      allocate (grid1_add_map1(max_links_map1), grid2_add_map1(max_links_map1), &
239                wts_map1(num_wts, max_links_map1))
240#ifdef TREAT_OVERLAY
241      allocate (grid1_add_repl1(grid1_size))
242#endif
243
244!-----------------------------------------------------------------------
245!
246!     allocate address and weight arrays for mapping 2 if necessary
247!   
248!-----------------------------------------------------------------------
249
250      if (num_maps > 1) then
251        allocate (grid1_add_map2(max_links_map2), &
252                  grid2_add_map2(max_links_map2), &
253                  wts_map2(num_wts, max_links_map2))
254      endif
255!-----------------------------------------------------------------------
256!
257!     initialize flag for routine store_link_cnsrv
258!   
259!-----------------------------------------------------------------------
260      first_call = .true.
261
262!-----------------------------------------------------------------------
263!
264      IF (nlogprt .GE. 2) THEN
265         WRITE (UNIT = nulou,FMT = *)' '
266         WRITE (UNIT = nulou,FMT = *)'Leaving routine init_remap_vars'
267         CALL OASIS_FLUSH_SCRIP(nulou)
268      ENDIF
269!
270      end subroutine init_remap_vars
271
272!***********************************************************************
273
274      subroutine thread_resize_remap_vars(self, increment)
275
276!-----------------------------------------------------------------------
277!
278!     this routine resizes remapping arrays by increasing(decreasing)
279!     the max_links by increment on a threaded structure
280!
281!-----------------------------------------------------------------------
282
283!-----------------------------------------------------------------------
284!
285!     input variables
286!
287!-----------------------------------------------------------------------
288
289      class(thread_remap), intent(inout) :: self
290      integer (kind=int_kind), intent(in) :: increment  ! the number of links to add(subtract) to arrays
291     
292!-----------------------------------------------------------------------
293!
294!     local variables
295!
296!-----------------------------------------------------------------------
297
298      integer (kind=int_kind) ::  ierr,    &  ! error flag
299                                  mxlinks   ! size of link arrays
300
301      integer (kind=int_kind), dimension(:), allocatable :: add1_tmp, &! temp array for resizing address arrays
302                                                            add2_tmp  ! temp array for resizing address arrays
303
304
305      real (kind=dbl_kind), dimension(:,:), allocatable :: wts_tmp   ! temp array for resizing weight arrays
306
307!***
308!*** allocate temporaries to hold original values
309!***
310
311      mxlinks = size(self%grid1_add)
312      allocate (add1_tmp(mxlinks), add2_tmp(mxlinks), wts_tmp(num_wts,mxlinks))
313
314      add1_tmp = self%grid1_add
315      add2_tmp = self%grid2_add
316      wts_tmp  = self%wts
317       
318!***
319!*** deallocate originals and increment max_links then
320!*** reallocate arrays at new size
321!***
322
323      deallocate (self%grid1_add, self%grid2_add, self%wts)
324      self%max_links = mxlinks + increment
325      allocate (self%grid1_add(self%max_links), self%grid2_add(self%max_links), self%wts(num_wts,self%max_links))
326
327!***
328!*** restore original values from temp arrays and
329!*** deallocate temps
330!***
331
332      mxlinks = min(mxlinks, self%max_links)
333      self%grid1_add(1:mxlinks) = add1_tmp (1:mxlinks)
334      self%grid2_add(1:mxlinks) = add2_tmp (1:mxlinks)
335      self%wts    (:,1:mxlinks) = wts_tmp(:,1:mxlinks)
336      deallocate(add1_tmp, add2_tmp, wts_tmp)
337
338      self%nb_resize = self%nb_resize + 1
339
340!
341!-----------------------------------------------------------------------
342!
343      end subroutine thread_resize_remap_vars
344!
345!-----------------------------------------------------------------------
346!-----------------------------------------------------------------------
347
348      subroutine resize_remap_vars(nmap, increment)
349
350!-----------------------------------------------------------------------
351!
352!     this routine resizes remapping arrays by increasing(decreasing)
353!     the max_links by increment
354!
355!-----------------------------------------------------------------------
356
357!-----------------------------------------------------------------------
358!
359!     input variables
360!
361!-----------------------------------------------------------------------
362
363      integer (kind=int_kind), intent(in) :: nmap,    &  ! identifies which mapping array to resize
364                                             increment  ! the number of links to add(subtract) to arrays
365
366!-----------------------------------------------------------------------
367!
368!     local variables
369!
370!-----------------------------------------------------------------------
371
372      integer (kind=int_kind) :: ierr,  &   ! error flag
373                                 mxlinks   ! size of link arrays
374
375      integer (kind=int_kind), dimension(:), allocatable :: add1_tmp, & ! temp array for resizing address arrays
376                                                            add2_tmp  ! temp array for resizing address arrays
377
378      real (kind=dbl_kind), dimension(:,:), allocatable ::  wts_tmp   ! temp array for resizing weight arrays
379!
380      IF (nlogprt .GE. 2) THEN
381         WRITE (UNIT = nulou,FMT = *)' '
382         WRITE (UNIT = nulou,FMT = *) 'Entering routine resize_remap_vars'
383         CALL OASIS_FLUSH_SCRIP(nulou)
384      ENDIF
385!
386!-----------------------------------------------------------------------
387!
388!     resize map 1 arrays if required.
389!
390!-----------------------------------------------------------------------
391
392      select case (nmap)
393      case(1)
394
395        !***
396        !*** allocate temporaries to hold original values
397        !***
398
399        mxlinks = size(grid1_add_map1)
400        allocate (add1_tmp(mxlinks), add2_tmp(mxlinks), wts_tmp(num_wts,mxlinks))
401
402        add1_tmp = grid1_add_map1
403        add2_tmp = grid2_add_map1
404        wts_tmp  = wts_map1
405       
406        !***
407        !*** deallocate originals and increment max_links then
408        !*** reallocate arrays at new size
409        !***
410
411        deallocate (grid1_add_map1, grid2_add_map1, wts_map1)
412        max_links_map1 = mxlinks + increment
413        allocate (grid1_add_map1(max_links_map1), grid2_add_map1(max_links_map1), wts_map1(num_wts,max_links_map1))
414
415        !***
416        !*** restore original values from temp arrays and
417        !*** deallocate temps
418        !***
419
420        mxlinks = min(mxlinks, max_links_map1)
421        grid1_add_map1(1:mxlinks) = add1_tmp (1:mxlinks)
422        grid2_add_map1(1:mxlinks) = add2_tmp (1:mxlinks)
423        wts_map1    (:,1:mxlinks) = wts_tmp(:,1:mxlinks)
424        deallocate(add1_tmp, add2_tmp, wts_tmp)
425
426!-----------------------------------------------------------------------
427!
428!     resize map 2 arrays if required.
429!
430!-----------------------------------------------------------------------
431
432      case(2)
433
434        !***
435        !*** allocate temporaries to hold original values
436        !***
437
438        mxlinks = size(grid1_add_map2)
439        allocate (add1_tmp(mxlinks), add2_tmp(mxlinks), wts_tmp(num_wts,mxlinks),stat=ierr)
440        if (ierr .ne. 0) THEN
441            WRITE(nulou,*) '   '
442            WRITE(nulou,*)'error allocating temps in resize: ',ierr
443            CALL OASIS_FLUSH_SCRIP(nulou)
444            stop
445        endif
446
447        add1_tmp = grid1_add_map2
448        add2_tmp = grid2_add_map2
449        wts_tmp  = wts_map2
450       
451        !***
452        !*** deallocate originals and increment max_links then
453        !*** reallocate arrays at new size
454        !***
455
456        deallocate (grid1_add_map2, grid2_add_map2, wts_map2)
457        max_links_map2 = mxlinks + increment
458        allocate (grid1_add_map2(max_links_map2), grid2_add_map2(max_links_map2), wts_map2(num_wts,max_links_map2),stat=ierr)
459        if (ierr .ne. 0) then
460            WRITE(nulou,*) '   '
461            WRITE(nulou,*)'error allocating new arrays in resize: ',ierr
462            CALL OASIS_FLUSH_SCRIP(nulou)
463            stop
464        endif
465
466
467        !***
468        !*** restore original values from temp arrays and
469        !*** deallocate temps
470        !***
471
472        mxlinks = min(mxlinks, max_links_map2)
473        grid1_add_map2(1:mxlinks) = add1_tmp (1:mxlinks)
474        grid2_add_map2(1:mxlinks) = add2_tmp (1:mxlinks)
475        wts_map2    (:,1:mxlinks) = wts_tmp(:,1:mxlinks)
476        deallocate(add1_tmp, add2_tmp, wts_tmp)
477
478      end select
479!
480      IF (nlogprt .GE. 2) THEN
481         WRITE (UNIT = nulou,FMT = *)' '
482         WRITE (UNIT = nulou,FMT = *) 'Leaving routine resize_remap_vars'
483         CALL OASIS_FLUSH_SCRIP(nulou)
484      ENDIF
485!
486!-----------------------------------------------------------------------
487
488      end subroutine resize_remap_vars
489!-----------------------------------------------------------------------
490!-----------------------------------------------------------------------
491      subroutine free_remap_vars
492!-----------------------------------------------------------------------
493!
494!     subroutine to deallocate allocated arrays
495!
496!-----------------------------------------------------------------------
497!
498      IF (nlogprt .GE. 2) THEN
499         WRITE (UNIT = nulou,FMT = *)' '
500         WRITE (UNIT = nulou,FMT = *)'Entering routine free_remap_vars'
501         CALL OASIS_FLUSH_SCRIP(nulou)
502      ENDIF
503!
504      deallocate (grid1_add_map1, grid2_add_map1, wts_map1)
505#ifdef TREAT_OVERLAY
506      deallocate (grid1_add_repl1)
507#endif
508
509      if (num_maps > 1) then
510        deallocate (grid1_add_map2, grid2_add_map2, wts_map2)
511      endif
512      if (map_type == map_type_conserv) then
513        first_call = .true.
514        first_conserv = .false.
515      endif
516!
517      IF (nlogprt .GE. 2) THEN
518         WRITE (UNIT = nulou,FMT = *)' '
519         WRITE (UNIT = nulou,FMT = *)'Leaving routine free_remap_vars'
520         CALL OASIS_FLUSH_SCRIP(nulou)
521      ENDIF
522!
523!-----------------------------------------------------------------------
524
525      end subroutine free_remap_vars
526
527!***********************************************************************
528
529      end module remap_vars
530
531!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Note: See TracBrowser for help on using the repository browser.