source: CPL/oasis3-mct_5.0/lib/scrip/src/fracnnei.F90 @ 6328

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

First import of oasis3-mct_5.0 (from oasis git server, branch OASIS3-MCT_5.0)
File size: 17.3 KB
Line 
1module fracnnei_mod
2
3contains
4
5   subroutine fracnnei()
6      !
7      !****
8      !               *****************************
9      !               * OASIS ROUTINE  -  LEVEL 4 *
10      !               * -------------     ------- *
11      !               *****************************
12      !
13      !**** *fracnnei* - SCRIP remapping
14      !
15      !     Purpose:
16      !     -------
17      !     Treatment of the tricky points in an interpolation
18      !
19      !     Interface:
20      !     ---------
21      !       *CALL*  *
22      !
23      !     Called from:
24      !     -----------
25      !     scriprmp
26      !
27      !     History:
28      !     -------
29      !       Version   Programmer     Date        Description
30      !       -------   ----------     ----        ----------- 
31      !       2.5       D.Declat       2002/08/20  adapted from S. Valcke ptmsq
32      !       3.0       S. Valcke      2002/10/30  test and corrections
33      !       4.0       A.Piacentini   2018/01/23  F90 and optimisation
34      !
35      ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
36      !* ---------------------------- Modules used ----------------------------
37      !
38      use kinds_mod     ! defines common data types
39      use constants     ! defines common constants
40      use grids         ! module containing grid information
41      use remap_vars    ! module containing remap information
42      use mod_oasis_flush
43!$ use omp_lib
44      !
45      !* ---------------------------- Implicit --------------------------------
46      !
47      implicit none
48      !
49      !* ---------------------------- Arguments -------------------------------
50      !
51      !
52      !* ---------------------------- Local Variables  ------------------------
53      !
54
55      logical (kind=log_kind) :: ll_debug = .false.
56
57      integer (kind=int_kind) ::  num_neigh, num_prev_links
58      logical (kind=log_kind), dimension(grid2_size) :: lnnei
59      integer (kind=int_kind) :: &
60         &    ib_dst,  &           ! INDEX loop for the distance grid
61         &    ib_src,  &           ! INDEX loop for the source grid
62         &    ib_links             ! INDEX loop for the links
63      integer (kind=int_kind) :: &
64         &    il_nneiadd, &        ! Nearest-neighbor address
65         &    il_lonkind, &
66         &    ib_bin,     &
67         &    min_add,    &
68         &    max_add
69      real (kind=dbl_kind) :: &
70         &    coslatd,   &          ! cosinus of the latitude for dst
71         &    sinlatd,   &          ! sinus of the latitude for dst
72         &    coslond,   &          ! cosinus of the longitude for dst
73         &    sinlond,   &          ! sinus of the longitude for dst
74         &    lonwestd,  &
75         &    loneastd,  &
76         &    latsouthd, &
77         &    latnorthd, &
78         &    distance, &
79         &    dist_min, &
80         &    d_dist
81
82      real (kind=dbl_kind), parameter :: dp_box     = 2.0*deg2rad
83      real (kind=dbl_kind), parameter :: dp_pthresh = 85.0*deg2rad
84
85      real (kind=dbl_kind), dimension(grid1_size) :: &
86         &    coslats,   &          ! cosinus of the latitude for src
87         &    sinlats,   &          ! sinus of the latitude for src
88         &    coslons,   &          ! cosinus of the longitude for src
89         &    sinlons               ! sinus of the longitude for src
90
91      integer (kind=int_kind) :: ib_vmm, il_add
92      integer (kind=int_kind), dimension(:), allocatable :: ila_dst
93      integer (kind=int_kind) :: il_envthreads, il_err, il_splitsize, ib_thread
94      integer (kind=int_kind) :: il_mythread
95      integer (kind=int_kind) :: il_nbthreads = 1
96      integer (kind=int_kind), dimension(:), allocatable :: ila_thr_sz
97      integer (kind=int_kind), dimension(:), allocatable :: ila_thr_mn
98      integer (kind=int_kind), dimension(:), allocatable :: ila_thr_mx
99
100      character (LEN=14) :: cl_envvar
101      !
102      !* ---------------------------- Poema verses ----------------------------
103      !
104      ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
105      !
106      !*    1. Initialization
107      !        --------------
108      !
109      if (nlogprt .ge. 2) then
110         write (UNIT = nulou,FMT = *) ' '
111         write (UNIT = nulou,FMT = *) ' '
112         write (UNIT = nulou,FMT = *) '   Entering ROUTINE fracnnei  -  Level 4'
113         write (UNIT = nulou,FMT = *) '           ****************     *******'
114         write (UNIT = nulou,FMT = *) ' '
115         write (UNIT = nulou,FMT = *) ' Treating the tricky points of the remapping'
116         write (UNIT = nulou,FMT = *) ' '
117         call OASIS_FLUSH_SCRIP(nulou)
118      endif
119      !
120      ! *----------------------------------------------------------------------
121      !
122      !*    2. Find Vmm points   V
123      !        ---------------  m m
124      !    A non-masked Valid target point is either with link or without link
125      !    If without link, find the non-masked nearest neighbours.
126      !
127      lnnei(:) = grid2_mask(:)
128      do ib_links = 1, num_links_map1
129         lnnei(grid2_add_map1(ib_links)) = .false.
130      end do
131
132      num_neigh = count(lnnei)
133
134      allocate(ila_dst(num_neigh))
135
136      ila_dst=pack([(ib_dst,ib_dst=1,grid2_size)],lnnei)
137
138      if (nlogprt .ge. 2) then
139         if (ll_debug) then
140            do ib_vmm = 1, num_neigh
141               write(nulou,*) '********* Will do FRACNNEI for point',ila_dst(ib_vmm)
142            end do
143         end if
144         write(nulou,*) '************ There are',num_neigh,'Vmm points in all'
145      end if
146
147      !
148      ! *----------------------------------------------------------------------
149      !
150      !*    2. Treating Vmm points   V
151      !        -------------------  m m
152      !     The target point is a non-masked Valid point while the source points
153      !         are all masked points. Use of the non-masked nearest neighbours.
154      !
155
156      !     Resize the storage
157
158      num_prev_links = num_links_map1
159      num_links_map1 = num_links_map1 + num_neigh
160
161      call resize_remap_vars(1, num_links_map1 - max_links_map1)
162
163      !* -- Precompute src grid spherical coordinates
164
165      coslats(:) = cos(grid1_center_lat(:))
166      coslons(:) = cos(grid1_center_lon(:))
167      sinlons(:) = sin(grid1_center_lon(:))
168      sinlats(:) = sin(grid1_center_lat(:))
169
170      call get_environment_variable(name='OASIS_OMP_NUM_THREADS', &
171         & value=cl_envvar, status=il_err)
172      if ( il_err .ne. 0) then
173         il_envthreads = 0
174      else
175         read(cl_envvar,*) il_envthreads
176      end if
177
178!$OMP PARALLEL NUM_THREADS(il_envthreads) DEFAULT(NONE) &
179!$OMP SHARED(il_envthreads)&
180!$OMP SHARED(sga_remap)&
181!$OMP SHARED(il_nbthreads) &
182!$OMP SHARED(ila_thr_sz,ila_thr_mn,ila_thr_mx) &
183!$OMP SHARED(num_neigh,ll_debug,ila_dst,num_prev_links,nulou) &
184!$OMP SHARED(num_wts) &
185!$OMP SHARED(grid1_add_map1,grid2_add_map1,wts_map1) &
186!$OMP SHARED(grid1_center_lat,grid1_center_lon) &
187!$OMP SHARED(grid2_center_lat,grid2_center_lon) &
188!$OMP SHARED(grid1_mask,grid1_size) &
189!$OMP SHARED(bin_lats,bin_addr1,num_srch_bins) &
190!$OMP SHARED(coslats,coslons,sinlons,sinlats) &
191!$OMP PRIVATE(nlogprt,ib_vmm,il_add,ib_dst,ib_src) &
192!$OMP PRIVATE(coslatd,sinlatd,coslond,sinlond) &
193!$OMP PRIVATE(lonwestd,loneastd,latsouthd,latnorthd) &
194!$OMP PRIVATE(dist_min,il_nneiadd,d_dist,distance) &
195!$OMP PRIVATE(il_lonkind,ib_bin,min_add,max_add) &
196!$OMP PRIVATE(ib_thread,il_splitsize)
197
198!$OMP SINGLE
199
200      il_nbthreads = 1
201!$    il_nbthreads = OMP_GET_NUM_THREADS ()
202
203      allocate(ila_thr_mn(il_nbthreads))
204      allocate(ila_thr_mx(il_nbthreads))
205
206      if (il_nbthreads .gt. 1) then
207
208         nlogprt = 0
209
210         allocate(ila_thr_sz(il_nbthreads))
211         ila_thr_sz(:) = floor(real(num_neigh)/il_nbthreads)
212         ila_thr_sz(1:num_neigh-sum(ila_thr_sz)) = ila_thr_sz(1:num_neigh-sum(ila_thr_sz)) + 1
213
214         ila_thr_mn(1) = 1
215         ila_thr_mx(1) = ila_thr_sz(1)
216
217         do ib_thread = 2, il_nbthreads
218            ila_thr_mn(ib_thread) = sum(ila_thr_sz(1:ib_thread-1)) + 1
219            ila_thr_mx(ib_thread) = sum(ila_thr_sz(1:ib_thread))
220         end do
221
222         allocate(sga_remap(il_nbthreads))
223
224         do ib_thread = 1, il_nbthreads
225            il_splitsize = ila_thr_sz(ib_thread)
226            sga_remap(ib_thread)%max_links = il_splitsize
227            sga_remap(ib_thread)%num_links = 0
228            allocate(sga_remap(ib_thread)%grid1_add(il_splitsize))
229            allocate(sga_remap(ib_thread)%grid2_add(il_splitsize))
230           !EM allocate(sga_remap(ib_thread)%wts(num_wts,il_splitsize))
231         end do
232
233         deallocate(ila_thr_sz)
234
235      else
236
237         ila_thr_mn(1) = 1
238         ila_thr_mx(1) = num_neigh
239
240      end if
241!$OMP END SINGLE
242
243      !* -- Find the nearest neighbours and store weights and address
244
245!$OMP DO SCHEDULE(STATIC,1)
246      thread_loop: do ib_thread = 1, il_nbthreads
247
248      grid_loop1: do ib_vmm = ila_thr_mn(ib_thread), ila_thr_mx(ib_thread)
249         ib_dst = ila_dst(ib_vmm)
250         il_add = num_prev_links+ib_vmm
251         if (il_nbthreads .eq. 1) then
252           grid2_add_map1(il_add) = ib_dst
253         else
254           sga_remap(ib_thread)%num_links = sga_remap(ib_thread)%num_links + 1
255           sga_remap(ib_thread)%grid2_add(sga_remap(ib_thread)%num_links) = ib_dst
256         endif
257
258         if (ll_debug) then
259            write(nulou,*) 'ib_dst for true=',ib_dst 
260            write(nulou,*) 'counter_Vmm =',ib_vmm
261            write(nulou,*) 'num_links+counter_Vmm =', il_add 
262            write(nulou,*) 'dst_addr =',grid2_add_map1(il_add)
263         endif
264
265         coslatd = cos(grid2_center_lat(ib_dst))
266         sinlatd = sin(grid2_center_lat(ib_dst))
267         coslond = cos(grid2_center_lon(ib_dst))
268         sinlond = sin(grid2_center_lon(ib_dst))
269
270         ! Draw a rough box around the target point (relaxed at poles)
271
272         if (grid2_center_lat(ib_dst) >= dp_pthresh) then
273
274            latsouthd  = grid2_center_lat(ib_dst) - dp_box
275            latnorthd  = pih
276            lonwestd   = zero
277            loneastd   = pi2
278            il_lonkind = 0
279
280         else if (grid2_center_lat(ib_dst) <= -dp_pthresh) then
281
282            latsouthd  = -pih
283            latnorthd  = grid2_center_lat(ib_dst) + dp_box
284            lonwestd   = zero
285            loneastd   = pi2
286            il_lonkind = 0
287
288         else
289
290            latsouthd = grid2_center_lat(ib_dst) - dp_box
291            latnorthd = grid2_center_lat(ib_dst) + dp_box
292            lonwestd  = grid2_center_lon(ib_dst) - dp_box/coslatd
293            loneastd  = grid2_center_lon(ib_dst) + dp_box/coslatd
294
295            if (lonwestd < zero) then
296               il_lonkind = -1
297            else if (loneastd > pi2) then
298               il_lonkind = 1
299            else
300               il_lonkind = 0
301            end if
302
303         end if
304
305         dist_min = bignum
306         il_nneiadd = 0
307
308         ! Restrict the loop by bins
309
310         min_add = grid1_size
311         max_add = 1
312         do ib_bin=1,num_srch_bins
313            if (bin_lats(2,ib_bin).ge.latsouthd .and.&
314              & bin_lats(1,ib_bin).le.latnorthd) then
315               min_add = min(min_add, bin_addr1(1,ib_bin))
316               max_add = max(max_add, bin_addr1(2,ib_bin))
317            end if
318         end do
319
320         do ib_src = min_add, max_add
321
322            ! Restrict the angular distance computation to a rough box
323
324            if (grid1_mask(ib_src)) then
325
326               if (grid1_center_lat(ib_src).le.latnorthd) then
327                  if (grid1_center_lat(ib_src).ge.latsouthd) then
328                     if (lf_loncheck(grid1_center_lon(ib_src),&
329                        &            lonwestd,loneastd,il_lonkind)) then
330                        d_dist = coslatd*coslats(ib_src) * &
331                           &  (coslond*coslons(ib_src) + &
332                           &   sinlond*sinlons(ib_src)) + &
333                           &  sinlatd*sinlats(ib_src)
334                        if (d_dist < -1.0d0) then
335                           d_dist = -1.0d0
336                        else if (d_dist > 1.0d0) then
337                           d_dist = 1.0d0
338                        end if
339                        distance = acos(d_dist)
340                        if (distance < dist_min) then
341                           il_nneiadd = ib_src
342                           dist_min = distance
343                        end if
344                     end if
345                  end if
346               end if
347            end if
348         end do
349
350         ! If the rough box was too small, search the entire grid
351
352         if (il_nneiadd == 0 ) then
353            dist_min = bignum
354            do ib_src = 1, grid1_size
355               if (grid1_mask(ib_src)) then
356                  d_dist = coslatd*coslats(ib_src) * &
357                     (coslond*coslons(ib_src) + &
358                     sinlond*sinlons(ib_src)) + &
359                     sinlatd*sinlats(ib_src)
360                  if (d_dist < -1.0d0) then
361                     d_dist = -1.0d0
362                  else if (d_dist > 1.0d0) then
363                     d_dist = 1.0d0
364                  end if
365                  distance = acos(d_dist)
366                  if (distance < dist_min) then
367                     il_nneiadd = ib_src
368                     dist_min = distance
369                  end if
370               end if
371            end do
372         end if
373
374         if (il_nbthreads .eq. 1) then
375           grid1_add_map1(il_add) = il_nneiadd
376           wts_map1(1,il_add) = 1.0
377           if (num_wts .gt. 1) then
378              wts_map1(2,il_add) = 0.0
379              wts_map1(3,il_add) = 0.0
380           endif
381         else
382           sga_remap(ib_thread)%grid1_add(sga_remap(ib_thread)%num_links) = il_nneiadd
383         endif
384         if (ll_debug.and.nlogprt .ge. 2) then
385            write(nulou,*) 'src_addr =',grid1_add_map1(il_add)
386            write(nulou,*) '*************** Nearest source neighbour is ',il_nneiadd         
387         end if
388      end do grid_loop1
389
390      end do thread_loop
391!$OMP END DO
392
393!$OMP END PARALLEL
394
395      if (il_nbthreads .gt. 1) then
396         sga_remap(1)%start_pos = num_prev_links + 1
397         il_splitsize = sga_remap(1)%num_links
398         do ib_thread = 2, il_nbthreads
399            il_splitsize = il_splitsize + sga_remap(ib_thread)%num_links
400            sga_remap(ib_thread)%start_pos = sga_remap(ib_thread-1)%start_pos + &
401               sga_remap(ib_thread-1)%num_links
402         end do
403
404         do ib_thread = 1, il_nbthreads
405            grid1_add_map1(sga_remap(ib_thread)%start_pos: &
406               sga_remap(ib_thread)%start_pos+             &
407               sga_remap(ib_thread)%num_links-1) =         &
408               sga_remap(ib_thread)%grid1_add
409            grid2_add_map1(sga_remap(ib_thread)%start_pos: &
410               sga_remap(ib_thread)%start_pos+             &
411               sga_remap(ib_thread)%num_links-1) =         &
412               sga_remap(ib_thread)%grid2_add
413            wts_map1     (1,sga_remap(ib_thread)%start_pos: &
414               sga_remap(ib_thread)%start_pos+            &
415               sga_remap(ib_thread)%num_links-1) = 1.0
416            wts_map1     (2,sga_remap(ib_thread)%start_pos: &
417               sga_remap(ib_thread)%start_pos+            &
418               sga_remap(ib_thread)%num_links-1) = 0.0
419            wts_map1     (3,sga_remap(ib_thread)%start_pos: &
420               sga_remap(ib_thread)%start_pos+            &
421               sga_remap(ib_thread)%num_links-1) = 0.0
422         end do
423
424         if (nlogprt.ge.2) then
425
426            do ib_thread = 1, il_nbthreads
427               if (sga_remap(ib_thread)%nb_resize.gt.0) then
428                  write(nulou,*) ' Number of thread_resize_remap_vars on thread ',&
429                     ib_thread, ' = ', sga_remap(ib_thread)%nb_resize
430               end if
431            end do
432
433         end if
434
435      end if
436
437      deallocate(ila_thr_mn)
438      deallocate(ila_thr_mx)
439      if (il_nbthreads .gt. 1) then
440         do ib_thread = 1, il_nbthreads
441            deallocate(sga_remap(ib_thread)%grid1_add)
442            deallocate(sga_remap(ib_thread)%grid2_add)
443            !EM deallocate(sga_remap(ib_thread)%wts)
444         end do
445         deallocate(sga_remap)
446      end if
447
448
449
450
451      !
452      deallocate(ila_dst)
453      !
454      ! *----------------------------------------------------------------------
455      !
456      if (nlogprt .ge. 2) then
457         write (UNIT = nulou,FMT = *) ' '
458         write (UNIT = nulou,FMT = *) '   Leaving ROUTINE fracnnei  -  Level 4'
459         write (UNIT = nulou,FMT = *) ' '
460         call OASIS_FLUSH_SCRIP(nulou)
461      end if
462
463   contains
464
465      logical (kind=log_kind) function lf_loncheck(rd_lon,rd_west,rd_east,id_kind)
466
467         real (kind=dbl_kind), intent(in) :: rd_lon
468         real (kind=dbl_kind), intent(in) :: rd_west,rd_east
469         integer (kind=int_kind), intent(in) :: id_kind
470         
471         select case(id_kind)
472
473         case(0)
474            lf_loncheck = (rd_lon>=rd_west) .and. (rd_lon<=rd_east)
475         case(-1)
476            lf_loncheck = (rd_lon-pi2>=rd_west) .or. (rd_lon<=rd_east)
477         case(1)
478            lf_loncheck = (rd_lon>=rd_west) .or. (rd_lon+pi2<=rd_east)
479         end select
480
481      end function lf_loncheck
482
483   end subroutine fracnnei
484
485   !***********************************************************************
486
487
488end module fracnnei_mod
Note: See TracBrowser for help on using the repository browser.