FV3 Bundle
sorted_index_nlm.F90
Go to the documentation of this file.
1 !***********************************************************************
2 !* GNU General Public License *
3 !* This file is a part of fvGFS. *
4 !* *
5 !* fvGFS is free software; you can redistribute it and/or modify it *
6 !* and are expected to follow the terms of the GNU General Public *
7 !* License as published by the Free Software Foundation; either *
8 !* version 2 of the License, or (at your option) any later version. *
9 !* *
10 !* fvGFS is distributed in the hope that it will be useful, but *
11 !* WITHOUT ANY WARRANTY; without even the implied warranty of *
12 !* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU *
13 !* General Public License for more details. *
14 !* *
15 !* For the full text of the GNU General Public License, *
16 !* write to: Free Software Foundation, Inc., *
17 !* 675 Mass Ave, Cambridge, MA 02139, USA. *
18 !* or see: http://www.gnu.org/licenses/gpl.html *
19 !***********************************************************************
20 !-*- F90 -*-
22  !---------------------------------------------------------------------
23  !<OVERVIEW>
24  ! sort cell corner indices in latlon space to ensure same order of
25  ! operations regardless of orientation in index space
26  !</OVERVIEW>
27  !
28  !<DESCRIPTION>
29  ! i/jinta are indices of b-grid locations needed for line integrals
30  ! around an a-grid cell including ghosting.
31  !
32  ! i/jintb are indices of a-grid locations needed for line integrals
33  ! around a b-grid cell, no ghosting.
34  !</DESCRIPTION>
35  !---------------------------------------------------------------------
36 
37  use fv_arrays_nlm_mod, only: r_grid
38 
39  implicit none
40  private
41  public :: sorted_inta, sorted_intb
42 
43 contains
44  !#####################################################################
45  ! <SUBROUTINE NAME="sorted_inta">
46  !
47  ! <DESCRIPTION>
48  ! Sort cell corner indices in latlon space based on grid locations
49  ! in index space. If not cubed_sphere assume orientations in index
50  ! and latlon space are identical.
51  !
52  ! i/jinta are indices of b-grid locations needed for line integrals
53  ! around an a-grid cell including ghosting.
54  !
55  ! i/jintb are indices of a-grid locations needed for line integrals
56  ! around a b-grid cell, no ghosting.
57  ! </DESCRIPTION>
58  !
59  subroutine sorted_inta(isd, ied, jsd, jed, cubed_sphere, bgrid, iinta, jinta)
60 
61  integer, intent(in) :: isd, ied, jsd, jed
62  real(kind=R_GRID), intent(in), dimension(isd:ied+1,jsd:jed+1,2) :: bgrid
63  logical, intent(in) :: cubed_sphere
64 
65  integer, intent(out), dimension(4,isd:ied,jsd:jed) :: iinta, jinta
66  !------------------------------------------------------------------!
67  ! local variables !
68  !------------------------------------------------------------------!
69  real, dimension(4) :: xsort, ysort
70  integer, dimension(4) :: isort, jsort
71  integer :: i, j
72  !------------------------------------------------------------------!
73  ! special treatment for cubed sphere !
74  !------------------------------------------------------------------!
75  if (cubed_sphere) then
76  !---------------------------------------------------------------!
77  ! get order of indices for line integral around a-grid cell !
78  !---------------------------------------------------------------!
79  do j=jsd,jed
80  do i=isd,ied
81  xsort(1)=bgrid(i ,j ,1); ysort(1)=bgrid(i ,j ,2); isort(1)=i ; jsort(1)=j
82  xsort(2)=bgrid(i ,j+1,1); ysort(2)=bgrid(i ,j+1,2); isort(2)=i ; jsort(2)=j+1
83  xsort(3)=bgrid(i+1,j+1,1); ysort(3)=bgrid(i+1,j+1,2); isort(3)=i+1; jsort(3)=j+1
84  xsort(4)=bgrid(i+1,j ,1); ysort(4)=bgrid(i+1,j ,2); isort(4)=i+1; jsort(4)=j
85  call sort_rectangle(iinta(1,i,j), jinta(1,i,j))
86  enddo
87  enddo
88  else
89  !---------------------------------------------------------------!
90  ! default behavior for other grids !
91  !---------------------------------------------------------------!
92  do j=jsd,jed
93  do i=isd,ied
94  iinta(i,j,1)=i ; jinta(i,j,1)=j
95  iinta(i,j,2)=i ; jinta(i,j,2)=j+1
96  iinta(i,j,3)=i+1; jinta(i,j,3)=j+1
97  iinta(i,j,4)=i+1; jinta(i,j,4)=j
98  enddo
99  enddo
100  endif
101 
102  contains
103  !------------------------------------------------------------------!
104  subroutine sort_rectangle(iind, jind)
105  integer, dimension(4), intent(inout) :: iind, jind
106  !----------------------------------------------------------------!
107  ! local variables !
108  !----------------------------------------------------------------!
109  real, dimension(4) :: xsorted, ysorted
110  integer, dimension(4) :: isorted, jsorted
111  integer :: l, ll, lll
112  !----------------------------------------------------------------!
113  ! sort in east west !
114  !----------------------------------------------------------------!
115  xsorted(:)=10.
116  ysorted(:)=10.
117  isorted(:)=0
118  jsorted(:)=0
119 
120  do l=1,4
121  do ll=1,4
122  if (xsort(l)<xsorted(ll)) then
123  do lll=3,ll,-1
124  xsorted(lll+1)=xsorted(lll)
125  ysorted(lll+1)=ysorted(lll)
126  isorted(lll+1)=isorted(lll)
127  jsorted(lll+1)=jsorted(lll)
128  enddo
129  xsorted(ll)=xsort(l)
130  ysorted(ll)=ysort(l)
131  isorted(ll)=isort(l)
132  jsorted(ll)=jsort(l)
133  exit
134  endif
135  enddo
136  enddo
137  !----------------------------------------------------------------!
138  ! sort in north south !
139  !----------------------------------------------------------------!
140  do l=1,4
141  xsort(l)=xsorted(l); ysort(l)=ysorted(l)
142  isort(l)=isorted(l); jsort(l)=jsorted(l)
143  enddo
144  xsorted(:)=10.
145  ysorted(:)=10.
146  isorted(:)=0
147  jsorted(:)=0
148 
149  do l=1,4
150  do ll=1,4
151  if (ysort(l)<ysorted(ll)) then
152  do lll=3,ll,-1
153  xsorted(lll+1)=xsorted(lll)
154  ysorted(lll+1)=ysorted(lll)
155  isorted(lll+1)=isorted(lll)
156  jsorted(lll+1)=jsorted(lll)
157  enddo
158  xsorted(ll)=xsort(l)
159  ysorted(ll)=ysort(l)
160  isorted(ll)=isort(l)
161  jsorted(ll)=jsort(l)
162  exit
163  endif
164  enddo
165  enddo
166  !----------------------------------------------------------------!
167  ! use first two grid point for start and orientation !
168  !----------------------------------------------------------------!
169  if ( isorted(1)==i .and. jsorted(1)==j ) then
170  if ( isorted(2)==i+1 .and. jsorted(2)==j+1 ) then
171  isorted(2)=isorted(3); jsorted(2)=jsorted(3)
172  endif
173  if ( isorted(2)==i .and. jsorted(2)==j+1 ) then
174  iind(1)=i ; jind(1)=j
175  iind(2)=i ; jind(2)=j+1
176  iind(3)=i+1; jind(3)=j+1
177  iind(4)=i+1; jind(4)=j
178  elseif ( isorted(2)==i+1 .and. jsorted(2)==j ) then
179  iind(1)=i ; jind(1)=j
180  iind(2)=i+1; jind(2)=j
181  iind(3)=i+1; jind(3)=j+1
182  iind(4)=i ; jind(4)=j+1
183  endif
184 
185  elseif ( isorted(1)==i .and. jsorted(1)==j+1 ) then
186  if ( isorted(2)==i+1 .and. jsorted(2)==j ) then
187  isorted(2)=isorted(3); jsorted(2)=jsorted(3)
188  endif
189  if ( isorted(2)==i+1 .and. jsorted(2)==j+1 ) then
190  iind(1)=i ; jind(1)=j+1
191  iind(2)=i+1; jind(2)=j+1
192  iind(3)=i+1; jind(3)=j
193  iind(4)=i ; jind(4)=j
194  elseif ( isorted(2)==i .and. jsorted(2)==j ) then
195  iind(1)=i ; jind(1)=j+1
196  iind(2)=i ; jind(2)=j
197  iind(3)=i+1; jind(3)=j
198  iind(4)=i+1; jind(4)=j+1
199  endif
200 
201  elseif ( isorted(1)==i+1 .and. jsorted(1)==j+1 ) then
202  if ( isorted(2)==i .and. jsorted(2)==j ) then
203  isorted(2)=isorted(3); jsorted(2)=jsorted(3)
204  endif
205  if ( isorted(2)==i+1 .and. jsorted(2)==j ) then
206  iind(1)=i+1; jind(1)=j+1
207  iind(2)=i+1; jind(2)=j
208  iind(3)=i ; jind(3)=j
209  iind(4)=i ; jind(4)=j+1
210  elseif ( isorted(2)==i .and. jsorted(2)==j+1 ) then
211  iind(1)=i+1; jind(1)=j+1
212  iind(2)=i ; jind(2)=j+1
213  iind(3)=i ; jind(3)=j
214  iind(4)=i+1; jind(4)=j
215  endif
216 
217  elseif ( isorted(1)==i+1 .and. jsorted(1)==j ) then
218  if ( isorted(2)==i .and. jsorted(2)==j+1 ) then
219  isorted(2)=isorted(3); jsorted(2)=jsorted(3)
220  endif
221  if ( isorted(2)==i .and. jsorted(2)==j ) then
222  iind(1)=i+1; jind(1)=j
223  iind(2)=i ; jind(2)=j
224  iind(3)=i ; jind(3)=j+1
225  iind(4)=i+1; jind(4)=j+1
226  elseif ( isorted(2)==i+1 .and. jsorted(2)==j+1 ) then
227  iind(1)=i+1; jind(1)=j
228  iind(2)=i+1; jind(2)=j+1
229  iind(3)=i ; jind(3)=j+1
230  iind(4)=i ; jind(4)=j
231  endif
232 
233  endif
234 
235  end subroutine sort_rectangle
236  !------------------------------------------------------------------!
237  end subroutine sorted_inta
238  ! </SUBROUTINE> NAME="sorted_inta"
239  !#####################################################################
240  ! <SUBROUTINE NAME="sorted_intb">
241  !
242  ! <DESCRIPTION>
243  ! Sort cell corner indices in latlon space based on grid locations
244  ! in index space. If not cubed_sphere assume orientations in index
245  ! and latlon space are identical.
246  !
247  ! i/jinta are indices of b-grid locations needed for line integrals
248  ! around an a-grid cell including ghosting.
249  !
250  ! i/jintb are indices of a-grid locations needed for line integrals
251  ! around a b-grid cell, no ghosting.
252  ! </DESCRIPTION>
253  !
254  subroutine sorted_intb(isd, ied, jsd, jed, is, ie, js, je, npx, npy, &
255  cubed_sphere, agrid, iintb, jintb)
257  integer, intent(in) :: isd, ied, jsd, jed, is, ie, js, je, npx, npy
258  real(kind=R_GRID), intent(in), dimension(isd:ied,jsd:jed,2) :: agrid
259  logical, intent(in) :: cubed_sphere
260 
261  integer, dimension(4,is:ie+1,js:je+1), intent(out) :: iintb, jintb
262  !------------------------------------------------------------------!
263  ! local variables !
264  !------------------------------------------------------------------!
265  real, dimension(4) :: xsort, ysort, xsorted, ysorted
266  integer, dimension(4) :: isort, jsort, isorted, jsorted
267  integer :: i, j, l, ll, lll
268  !------------------------------------------------------------------!
269  ! special treatment for cubed sphere !
270  !------------------------------------------------------------------!
271  if (cubed_sphere) then
272  !---------------------------------------------------------------!
273  ! get order of indices for line integral around b-grid cell !
274  !---------------------------------------------------------------!
275  do j=js,je+1
276  do i=is,ie+1
277  xsort(1)=agrid(i ,j ,1); ysort(1)=agrid(i ,j ,2); isort(1)=i ; jsort(1)=j
278  xsort(2)=agrid(i ,j-1,1); ysort(2)=agrid(i ,j-1,2); isort(2)=i ; jsort(2)=j-1
279  xsort(3)=agrid(i-1,j-1,1); ysort(3)=agrid(i-1,j-1,2); isort(3)=i-1; jsort(3)=j-1
280  xsort(4)=agrid(i-1,j ,1); ysort(4)=agrid(i-1,j ,2); isort(4)=i-1; jsort(4)=j
281  call sort_rectangle(iintb(1,i,j), jintb(1,i,j))
282  enddo
283  enddo
284  !---------------------------------------------------------------!
285  ! take care of corner points !
286  !---------------------------------------------------------------!
287  if ( (is==1) .and. (js==1) ) then
288  i=1
289  j=1
290  xsort(1)=agrid(i ,j ,1); ysort(1)=agrid(i ,j ,2); isort(1)=i ; jsort(1)=j
291  xsort(2)=agrid(i ,j-1,1); ysort(2)=agrid(i ,j-1,2); isort(2)=i ; jsort(2)=j-1
292  xsort(3)=agrid(i-1,j ,1); ysort(3)=agrid(i-1,j ,2); isort(3)=i-1; jsort(3)=j
293  call sort_triangle()
294  iintb(4,i,j)=i-1; jintb(4,i,j)=j-1
295  endif
296 
297  if ( (ie+1==npx) .and. (js==1) ) then
298  i=npx
299  j=1
300  xsort(1)=agrid(i ,j ,1); ysort(1)=agrid(i ,j ,2); isort(1)=i ; jsort(1)=j
301  xsort(2)=agrid(i-1,j ,1); ysort(2)=agrid(i-1,j ,2); isort(2)=i-1; jsort(2)=j
302  xsort(3)=agrid(i-1,j-1,1); ysort(3)=agrid(i-1,j-1,2); isort(3)=i-1; jsort(3)=j-1
303  call sort_triangle()
304  iintb(4,i,j)=i; jintb(4,i,j)=j-1
305  endif
306 
307  if ( (ie+1==npx) .and. (je+1==npy) ) then
308  i=npx
309  j=npy
310  xsort(1)=agrid(i-1,j-1,1); ysort(1)=agrid(i-1,j-1,2); isort(1)=i-1; jsort(1)=j-1
311  xsort(2)=agrid(i ,j-1,1); ysort(2)=agrid(i ,j-1,2); isort(2)=i ; jsort(2)=j-1
312  xsort(3)=agrid(i-1,j ,1); ysort(3)=agrid(i-1,j ,2); isort(3)=i-1; jsort(3)=j
313  call sort_triangle()
314  iintb(4,i,j)=i; jintb(4,i,j)=j
315  endif
316 
317  if ( (is==1) .and. (je+1==npy) ) then
318  i=1
319  j=npy
320  xsort(1)=agrid(i ,j ,1); ysort(1)=agrid(i ,j ,2); isort(1)=i ; jsort(1)=j
321  xsort(2)=agrid(i-1,j-1,1); ysort(2)=agrid(i-1,j-1,2); isort(2)=i-1; jsort(2)=j-1
322  xsort(3)=agrid(i ,j-1,1); ysort(3)=agrid(i ,j-1,2); isort(3)=i ; jsort(3)=j-1
323  call sort_triangle()
324  iintb(4,i,j)=i-1; jintb(4,i,j)=j
325  endif
326  else
327  !---------------------------------------------------------------!
328  ! default behavior for other grids !
329  !---------------------------------------------------------------!
330  do j=js,je+1
331  do i=is,ie+1
332  iintb(1,i,j)=i ; jintb(1,i,j)=j
333  iintb(2,i,j)=i ; jintb(2,i,j)=j-1
334  iintb(3,i,j)=i-1; jintb(3,i,j)=j-1
335  iintb(4,i,j)=i-1; jintb(4,i,j)=j
336  enddo
337  enddo
338  endif
339 
340  contains
341  !------------------------------------------------------------------!
342  subroutine sort_rectangle(iind, jind)
343 
344  integer, dimension(4), intent(inout) :: iind, jind
345  !----------------------------------------------------------------!
346  ! local variables !
347  !----------------------------------------------------------------!
348  real, dimension(4) :: xsorted, ysorted
349  integer, dimension(4) :: isorted, jsorted
350  !----------------------------------------------------------------!
351  ! sort in east west !
352  !----------------------------------------------------------------!
353  xsorted(:)=10.
354  ysorted(:)=10.
355  isorted(:)=0
356  jsorted(:)=0
357 
358  do l=1,4
359  do ll=1,4
360  if (xsort(l)<xsorted(ll)) then
361  do lll=3,ll,-1
362  xsorted(lll+1)=xsorted(lll)
363  ysorted(lll+1)=ysorted(lll)
364  isorted(lll+1)=isorted(lll)
365  jsorted(lll+1)=jsorted(lll)
366  enddo
367  xsorted(ll)=xsort(l)
368  ysorted(ll)=ysort(l)
369  isorted(ll)=isort(l)
370  jsorted(ll)=jsort(l)
371  exit
372  endif
373  enddo
374  enddo
375  !----------------------------------------------------------------!
376  ! sort in north south !
377  !----------------------------------------------------------------!
378  do l=1,4
379  xsort(l)=xsorted(l); ysort(l)=ysorted(l)
380  isort(l)=isorted(l); jsort(l)=jsorted(l)
381  enddo
382  xsorted(:)=10.
383  ysorted(:)=10.
384  isorted(:)=0
385  jsorted(:)=0
386 
387  do l=1,4
388  do ll=1,4
389  if (ysort(l)<ysorted(ll)) then
390  do lll=3,ll,-1
391  xsorted(lll+1)=xsorted(lll)
392  ysorted(lll+1)=ysorted(lll)
393  isorted(lll+1)=isorted(lll)
394  jsorted(lll+1)=jsorted(lll)
395  enddo
396  xsorted(ll)=xsort(l)
397  ysorted(ll)=ysort(l)
398  isorted(ll)=isort(l)
399  jsorted(ll)=jsort(l)
400  exit
401  endif
402  enddo
403  enddo
404  !----------------------------------------------------------------!
405  ! use first two grid point for start and orientation !
406  !----------------------------------------------------------------!
407  if ( isorted(1)==i .and. jsorted(1)==j ) then
408  if ( isorted(2)==i-1 .and. jsorted(2)==j-1 ) then
409  isorted(2)=isorted(3); jsorted(2)=jsorted(3)
410  endif
411  if ( isorted(2)==i .and. jsorted(2)==j-1 ) then
412  iind(1)=i ; jind(1)=j
413  iind(2)=i ; jind(2)=j-1
414  iind(3)=i-1; jind(3)=j-1
415  iind(4)=i-1; jind(4)=j
416  elseif ( isorted(2)==i-1 .and. jsorted(2)==j ) then
417  iind(1)=i ; jind(1)=j
418  iind(2)=i-1; jind(2)=j
419  iind(3)=i-1; jind(3)=j-1
420  iind(4)=i ; jind(4)=j-1
421  endif
422 
423  elseif ( isorted(1)==i .and. jsorted(1)==j-1 ) then
424  if ( isorted(2)==i-1 .and. jsorted(2)==j ) then
425  isorted(2)=isorted(3); jsorted(2)=jsorted(3)
426  endif
427  if ( isorted(2)==i-1 .and. jsorted(2)==j-1 ) then
428  iind(1)=i ; jind(1)=j-1
429  iind(2)=i-1; jind(2)=j-1
430  iind(3)=i-1; jind(3)=j
431  iind(4)=i ; jind(4)=j
432  elseif ( isorted(2)==i .and. jsorted(2)==j ) then
433  iind(1)=i ; jind(1)=j-1
434  iind(2)=i ; jind(2)=j
435  iind(3)=i-1; jind(3)=j
436  iind(4)=i-1; jind(4)=j-1
437  endif
438 
439  elseif ( isorted(1)==i-1 .and. jsorted(1)==j-1 ) then
440  if ( isorted(2)==i .and. jsorted(2)==j ) then
441  isorted(2)=isorted(3); jsorted(2)=jsorted(3)
442  endif
443  if ( isorted(2)==i-1 .and. jsorted(2)==j ) then
444  iind(1)=i-1; jind(1)=j-1
445  iind(2)=i-1; jind(2)=j
446  iind(3)=i ; jind(3)=j
447  iind(4)=i ; jind(4)=j-1
448  elseif ( isorted(2)==i .and. jsorted(2)==j-1 ) then
449  iind(1)=i-1; jind(1)=j-1
450  iind(2)=i ; jind(2)=j-1
451  iind(3)=i ; jind(3)=j
452  iind(4)=i-1; jind(4)=j
453  endif
454 
455  elseif ( isorted(1)==i-1 .and. jsorted(1)==j ) then
456  if ( isorted(2)==i .and. jsorted(2)==j-1 ) then
457  isorted(2)=isorted(3); jsorted(2)=jsorted(3)
458  endif
459  if ( isorted(2)==i .and. jsorted(2)==j ) then
460  iind(1)=i-1; jind(1)=j
461  iind(2)=i ; jind(2)=j
462  iind(3)=i ; jind(3)=j-1
463  iind(4)=i-1; jind(4)=j-1
464  elseif ( isorted(2)==i-1 .and. jsorted(2)==j-1 ) then
465  iind(1)=i-1; jind(1)=j
466  iind(2)=i-1; jind(2)=j-1
467  iind(3)=i ; jind(3)=j-1
468  iind(4)=i ; jind(4)=j
469  endif
470 
471  endif
472 
473  end subroutine sort_rectangle
474  !------------------------------------------------------------------!
475  subroutine sort_triangle()
477  xsorted(1:3)=10.
478  ysorted(1:3)=10.
479  isorted(1:3)=0
480  jsorted(1:3)=0
481  !----------------------------------------------------------------!
482  ! sort in east west !
483  !----------------------------------------------------------------!
484  do l=1,3
485  do ll=1,3
486  if (xsort(l)<xsorted(ll)) then
487  do lll=2,ll,-1
488  xsorted(lll+1)=xsorted(lll)
489  ysorted(lll+1)=ysorted(lll)
490  isorted(lll+1)=isorted(lll)
491  jsorted(lll+1)=jsorted(lll)
492  enddo
493  xsorted(ll)=xsort(l)
494  ysorted(ll)=ysort(l)
495  isorted(ll)=isort(l)
496  jsorted(ll)=jsort(l)
497  exit
498  endif
499  enddo
500  enddo
501  !----------------------------------------------------------------!
502  ! sort in north south !
503  !----------------------------------------------------------------!
504  do l=1,3
505  xsort(l)=xsorted(l); ysort(l)=ysorted(l)
506  isort(l)=isorted(l); jsort(l)=jsorted(l)
507  enddo
508  xsorted(1:3)=10.
509  ysorted(1:3)=10.
510  isorted(1:3)=0
511  jsorted(1:3)=0
512 
513  do l=1,3
514  do ll=1,3
515  if (ysort(l)<ysorted(ll)) then
516  do lll=2,ll,-1
517  xsorted(lll+1)=xsorted(lll)
518  ysorted(lll+1)=ysorted(lll)
519  isorted(lll+1)=isorted(lll)
520  jsorted(lll+1)=jsorted(lll)
521  enddo
522  xsorted(ll)=xsort(l)
523  ysorted(ll)=ysort(l)
524  isorted(ll)=isort(l)
525  jsorted(ll)=jsort(l)
526  exit
527  endif
528  enddo
529  enddo
530  !----------------------------------------------------------------!
531  ! use first two grid point for start and orientation !
532  !----------------------------------------------------------------!
533  iintb(1,i,j)=isorted(1) ; jintb(1,i,j)=jsorted(1)
534  iintb(2,i,j)=isorted(2) ; jintb(2,i,j)=jsorted(2)
535  iintb(3,i,j)=isorted(3) ; jintb(3,i,j)=jsorted(3)
536 
537  end subroutine sort_triangle
538  !------------------------------------------------------------------!
539  end subroutine sorted_intb
540  ! </SUBROUTINE> NAME="sorted_intb"
541  !#####################################################################
542 end module sorted_index_nlm_mod
integer, parameter, public r_grid
subroutine, public sorted_intb(isd, ied, jsd, jed, is, ie, js, je, npx, npy, cubed_sphere, agrid, iintb, jintb)
subroutine sort_triangle()
subroutine sort_rectangle(iind, jind)
subroutine, public sorted_inta(isd, ied, jsd, jed, cubed_sphere, bgrid, iinta, jinta)