source: LMDZ6/trunk/libf/phylmd/cosp/m_mrgrnk.f90 @ 5300

Last change on this file since 5300 was 5268, checked in by abarral, 5 weeks ago

.f90 <-> .F90 depending on cpp key use
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 11.2 KB
Line 
1Module m_mrgrnk
2IMPLICIT NONE
3Integer, Parameter :: kdp = selected_real_kind(15)
4public :: mrgrnk
5private :: kdp
6private :: I_mrgrnk, D_mrgrnk
7interface mrgrnk
8  module procedure D_mrgrnk, I_mrgrnk
9end interface mrgrnk
10contains
11
12Subroutine D_mrgrnk (XDONT, IRNGT)
13! __________________________________________________________
14!   MRGRNK = Merge-sort ranking of an array
15!   For performance reasons, the first 2 passes are taken
16!   out of the standard loop, and use dedicated coding.
17! __________________________________________________________
18! __________________________________________________________
19      Real (kind=kdp), Dimension (:), Intent (In) :: XDONT
20      Integer, Dimension (:), Intent (Out) :: IRNGT
21! __________________________________________________________
22      Real (kind=kdp) :: XVALA, XVALB
23!
24      Integer, Dimension (SIZE(IRNGT)) :: JWRKT
25      Integer :: LMTNA, LMTNC, IRNG1, IRNG2
26      Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
27!
28      NVAL = Min (SIZE(XDONT), SIZE(IRNGT))
29      Select Case (NVAL)
30      Case (:0)
31         Return
32      Case (1)
33         IRNGT (1) = 1
34         Return
35      Case Default
36         Continue
37      End Select
38!
39!  Fill-in the index array, creating ordered couples
40!
41      Do IIND = 2, NVAL, 2
42         If (XDONT(IIND-1) <= XDONT(IIND)) Then
43            IRNGT (IIND-1) = IIND - 1
44            IRNGT (IIND) = IIND
45         Else
46            IRNGT (IIND-1) = IIND
47            IRNGT (IIND) = IIND - 1
48         End If
49      End Do
50      If (Modulo(NVAL, 2) /= 0) Then
51         IRNGT (NVAL) = NVAL
52      End If
53!
54!  We will now have ordered subsets A - B - A - B - ...
55!  and merge A and B couples into     C   -   C   - ...
56!
57      LMTNA = 2
58      LMTNC = 4
59!
60!  First iteration. The length of the ordered subsets goes from 2 to 4
61!
62      Do
63         If (NVAL <= 2) Exit
64!
65!   Loop on merges of A and B into C
66!
67         Do IWRKD = 0, NVAL - 1, 4
68            If ((IWRKD+4) > NVAL) Then
69               If ((IWRKD+2) >= NVAL) Exit
70!
71!   1 2 3
72!
73               If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
74!
75!   1 3 2
76!
77               If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
78                  IRNG2 = IRNGT (IWRKD+2)
79                  IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
80                  IRNGT (IWRKD+3) = IRNG2
81!
82!   3 1 2
83!
84               Else
85                  IRNG1 = IRNGT (IWRKD+1)
86                  IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
87                  IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
88                  IRNGT (IWRKD+2) = IRNG1
89               End If
90               Exit
91            End If
92!
93!   1 2 3 4
94!
95            If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
96!
97!   1 3 x x
98!
99            If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
100               IRNG2 = IRNGT (IWRKD+2)
101               IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
102               If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
103!   1 3 2 4
104                  IRNGT (IWRKD+3) = IRNG2
105               Else
106!   1 3 4 2
107                  IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
108                  IRNGT (IWRKD+4) = IRNG2
109               End If
110!
111!   3 x x x
112!
113            Else
114               IRNG1 = IRNGT (IWRKD+1)
115               IRNG2 = IRNGT (IWRKD+2)
116               IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
117               If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
118                  IRNGT (IWRKD+2) = IRNG1
119                  If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
120!   3 1 2 4
121                     IRNGT (IWRKD+3) = IRNG2
122                  Else
123!   3 1 4 2
124                     IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
125                     IRNGT (IWRKD+4) = IRNG2
126                  End If
127               Else
128!   3 4 1 2
129                  IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
130                  IRNGT (IWRKD+3) = IRNG1
131                  IRNGT (IWRKD+4) = IRNG2
132               End If
133            End If
134         End Do
135!
136!  The Cs become As and Bs
137!
138         LMTNA = 4
139         Exit
140      End Do
141!
142!  Iteration loop. Each time, the length of the ordered subsets
143!  is doubled.
144!
145      Do
146         If (LMTNA >= NVAL) Exit
147         IWRKF = 0
148         LMTNC = 2 * LMTNC
149!
150!   Loop on merges of A and B into C
151!
152         Do
153            IWRK = IWRKF
154            IWRKD = IWRKF + 1
155            JINDA = IWRKF + LMTNA
156            IWRKF = IWRKF + LMTNC
157            If (IWRKF >= NVAL) Then
158               If (JINDA >= NVAL) Exit
159               IWRKF = NVAL
160            End If
161            IINDA = 1
162            IINDB = JINDA + 1
163!
164!   Shortcut for the case when the max of A is smaller
165!   than the min of B. This line may be activated when the
166!   initial set is already close to sorted.
167!
168!          IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
169!
170!  One steps in the C subset, that we build in the final rank array
171!
172!  Make a copy of the rank array for the merge iteration
173!
174            JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
175!
176            XVALA = XDONT (JWRKT(IINDA))
177            XVALB = XDONT (IRNGT(IINDB))
178!
179            Do
180               IWRK = IWRK + 1
181!
182!  We still have unprocessed values in both A and B
183!
184               If (XVALA > XVALB) Then
185                  IRNGT (IWRK) = IRNGT (IINDB)
186                  IINDB = IINDB + 1
187                  If (IINDB > IWRKF) Then
188!  Only A still with unprocessed values
189                     IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
190                     Exit
191                  End If
192                  XVALB = XDONT (IRNGT(IINDB))
193               Else
194                  IRNGT (IWRK) = JWRKT (IINDA)
195                  IINDA = IINDA + 1
196                  If (IINDA > LMTNA) Exit! Only B still with unprocessed values
197                  XVALA = XDONT (JWRKT(IINDA))
198               End If
199!
200            End Do
201         End Do
202!
203!  The Cs become As and Bs
204!
205         LMTNA = 2 * LMTNA
206      End Do
207!
208      Return
209!
210End Subroutine D_mrgrnk
211
212Subroutine I_mrgrnk (XDONT, IRNGT)
213! __________________________________________________________
214!   MRGRNK = Merge-sort ranking of an array
215!   For performance reasons, the first 2 passes are taken
216!   out of the standard loop, and use dedicated coding.
217! __________________________________________________________
218! __________________________________________________________
219      Integer, Dimension (:), Intent (In)  :: XDONT
220      Integer, Dimension (:), Intent (Out) :: IRNGT
221! __________________________________________________________
222      Integer :: XVALA, XVALB
223!
224      Integer, Dimension (SIZE(IRNGT)) :: JWRKT
225      Integer :: LMTNA, LMTNC, IRNG1, IRNG2
226      Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
227!
228      NVAL = Min (SIZE(XDONT), SIZE(IRNGT))
229      Select Case (NVAL)
230      Case (:0)
231         Return
232      Case (1)
233         IRNGT (1) = 1
234         Return
235      Case Default
236         Continue
237      End Select
238!
239!  Fill-in the index array, creating ordered couples
240!
241      Do IIND = 2, NVAL, 2
242         If (XDONT(IIND-1) <= XDONT(IIND)) Then
243            IRNGT (IIND-1) = IIND - 1
244            IRNGT (IIND) = IIND
245         Else
246            IRNGT (IIND-1) = IIND
247            IRNGT (IIND) = IIND - 1
248         End If
249      End Do
250      If (Modulo(NVAL, 2) /= 0) Then
251         IRNGT (NVAL) = NVAL
252      End If
253!
254!  We will now have ordered subsets A - B - A - B - ...
255!  and merge A and B couples into     C   -   C   - ...
256!
257      LMTNA = 2
258      LMTNC = 4
259!
260!  First iteration. The length of the ordered subsets goes from 2 to 4
261!
262      Do
263         If (NVAL <= 2) Exit
264!
265!   Loop on merges of A and B into C
266!
267         Do IWRKD = 0, NVAL - 1, 4
268            If ((IWRKD+4) > NVAL) Then
269               If ((IWRKD+2) >= NVAL) Exit
270!
271!   1 2 3
272!
273               If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
274!
275!   1 3 2
276!
277               If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
278                  IRNG2 = IRNGT (IWRKD+2)
279                  IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
280                  IRNGT (IWRKD+3) = IRNG2
281!
282!   3 1 2
283!
284               Else
285                  IRNG1 = IRNGT (IWRKD+1)
286                  IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
287                  IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
288                  IRNGT (IWRKD+2) = IRNG1
289               End If
290               Exit
291            End If
292!
293!   1 2 3 4
294!
295            If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
296!
297!   1 3 x x
298!
299            If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
300               IRNG2 = IRNGT (IWRKD+2)
301               IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
302               If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
303!   1 3 2 4
304                  IRNGT (IWRKD+3) = IRNG2
305               Else
306!   1 3 4 2
307                  IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
308                  IRNGT (IWRKD+4) = IRNG2
309               End If
310!
311!   3 x x x
312!
313            Else
314               IRNG1 = IRNGT (IWRKD+1)
315               IRNG2 = IRNGT (IWRKD+2)
316               IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
317               If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
318                  IRNGT (IWRKD+2) = IRNG1
319                  If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
320!   3 1 2 4
321                     IRNGT (IWRKD+3) = IRNG2
322                  Else
323!   3 1 4 2
324                     IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
325                     IRNGT (IWRKD+4) = IRNG2
326                  End If
327               Else
328!   3 4 1 2
329                  IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
330                  IRNGT (IWRKD+3) = IRNG1
331                  IRNGT (IWRKD+4) = IRNG2
332               End If
333            End If
334         End Do
335!
336!  The Cs become As and Bs
337!
338         LMTNA = 4
339         Exit
340      End Do
341!
342!  Iteration loop. Each time, the length of the ordered subsets
343!  is doubled.
344!
345      Do
346         If (LMTNA >= NVAL) Exit
347         IWRKF = 0
348         LMTNC = 2 * LMTNC
349!
350!   Loop on merges of A and B into C
351!
352         Do
353            IWRK = IWRKF
354            IWRKD = IWRKF + 1
355            JINDA = IWRKF + LMTNA
356            IWRKF = IWRKF + LMTNC
357            If (IWRKF >= NVAL) Then
358               If (JINDA >= NVAL) Exit
359               IWRKF = NVAL
360            End If
361            IINDA = 1
362            IINDB = JINDA + 1
363!
364!   Shortcut for the case when the max of A is smaller
365!   than the min of B. This line may be activated when the
366!   initial set is already close to sorted.
367!
368!          IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
369!
370!  One steps in the C subset, that we build in the final rank array
371!
372!  Make a copy of the rank array for the merge iteration
373!
374            JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
375!
376            XVALA = XDONT (JWRKT(IINDA))
377            XVALB = XDONT (IRNGT(IINDB))
378!
379            Do
380               IWRK = IWRK + 1
381!
382!  We still have unprocessed values in both A and B
383!
384               If (XVALA > XVALB) Then
385                  IRNGT (IWRK) = IRNGT (IINDB)
386                  IINDB = IINDB + 1
387                  If (IINDB > IWRKF) Then
388!  Only A still with unprocessed values
389                     IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
390                     Exit
391                  End If
392                  XVALB = XDONT (IRNGT(IINDB))
393               Else
394                  IRNGT (IWRK) = JWRKT (IINDA)
395                  IINDA = IINDA + 1
396                  If (IINDA > LMTNA) Exit! Only B still with unprocessed values
397                  XVALA = XDONT (JWRKT(IINDA))
398               End If
399!
400            End Do
401         End Do
402!
403!  The Cs become As and Bs
404!
405         LMTNA = 2 * LMTNA
406      End Do
407!
408      Return
409!
410End Subroutine I_mrgrnk
411end module m_mrgrnk
Note: See TracBrowser for help on using the repository browser.