source: LMDZ5/branches/AI-cosp/libf/phylmd/cosp/mrgrnk.F90.new @ 5440

Last change on this file since 5440 was 2432, checked in by idelkadi, 9 years ago

Correction : rajout de fichiers manquants pour les simulateurs Modis, Cloudsat, Misr et Isccp
File size: 16.8 KB
Line 
1! $Revision: 23 $, $Date: 2011-03-31 15:41:37 +0200 (jeu. 31 mars 2011) $
2! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/quickbeam/mrgrnk.f90 $
3Module m_mrgrnk
4Integer, Parameter :: kdp = selected_real_kind(15)
5public :: mrgrnk
6private :: kdp
7private :: R_mrgrnk, I_mrgrnk, D_mrgrnk
8interface mrgrnk
9  module procedure D_mrgrnk, R_mrgrnk, I_mrgrnk
10end interface mrgrnk
11contains
12
13Subroutine D_mrgrnk (XDONT, IRNGT)
14! __________________________________________________________
15!   MRGRNK = Merge-sort ranking of an array
16!   For performance reasons, the first 2 passes are taken
17!   out of the standard loop, and use dedicated coding.
18! __________________________________________________________
19! __________________________________________________________
20      Real (kind=kdp), Dimension (:), Intent (In) :: XDONT
21      Integer*4, Dimension (:), Intent (Out) :: IRNGT
22! __________________________________________________________
23      Real (kind=kdp) :: XVALA, XVALB
24!
25      Integer, Dimension (SIZE(IRNGT)) :: JWRKT
26      Integer :: LMTNA, LMTNC, IRNG1, IRNG2
27      Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
28!
29      NVAL = Min (SIZE(XDONT), SIZE(IRNGT))
30      Select Case (NVAL)
31      Case (:0)
32         Return
33      Case (1)
34         IRNGT (1) = 1
35         Return
36      Case Default
37         Continue
38      End Select
39!
40!  Fill-in the index array, creating ordered couples
41!
42      Do IIND = 2, NVAL, 2
43         If (XDONT(IIND-1) <= XDONT(IIND)) Then
44            IRNGT (IIND-1) = IIND - 1
45            IRNGT (IIND) = IIND
46         Else
47            IRNGT (IIND-1) = IIND
48            IRNGT (IIND) = IIND - 1
49         End If
50      End Do
51      If (Modulo(NVAL, 2) /= 0) Then
52         IRNGT (NVAL) = NVAL
53      End If
54!
55!  We will now have ordered subsets A - B - A - B - ...
56!  and merge A and B couples into     C   -   C   - ...
57!
58      LMTNA = 2
59      LMTNC = 4
60!
61!  First iteration. The length of the ordered subsets goes from 2 to 4
62!
63      Do
64         If (NVAL <= 2) Exit
65!
66!   Loop on merges of A and B into C
67!
68         Do IWRKD = 0, NVAL - 1, 4
69            If ((IWRKD+4) > NVAL) Then
70               If ((IWRKD+2) >= NVAL) Exit
71!
72!   1 2 3
73!
74               If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
75!
76!   1 3 2
77!
78               If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
79                  IRNG2 = IRNGT (IWRKD+2)
80                  IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
81                  IRNGT (IWRKD+3) = IRNG2
82!
83!   3 1 2
84!
85               Else
86                  IRNG1 = IRNGT (IWRKD+1)
87                  IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
88                  IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
89                  IRNGT (IWRKD+2) = IRNG1
90               End If
91               Exit
92            End If
93!
94!   1 2 3 4
95!
96            If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
97!
98!   1 3 x x
99!
100            If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
101               IRNG2 = IRNGT (IWRKD+2)
102               IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
103               If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
104!   1 3 2 4
105                  IRNGT (IWRKD+3) = IRNG2
106               Else
107!   1 3 4 2
108                  IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
109                  IRNGT (IWRKD+4) = IRNG2
110               End If
111!
112!   3 x x x
113!
114            Else
115               IRNG1 = IRNGT (IWRKD+1)
116               IRNG2 = IRNGT (IWRKD+2)
117               IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
118               If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
119                  IRNGT (IWRKD+2) = IRNG1
120                  If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
121!   3 1 2 4
122                     IRNGT (IWRKD+3) = IRNG2
123                  Else
124!   3 1 4 2
125                     IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
126                     IRNGT (IWRKD+4) = IRNG2
127                  End If
128               Else
129!   3 4 1 2
130                  IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
131                  IRNGT (IWRKD+3) = IRNG1
132                  IRNGT (IWRKD+4) = IRNG2
133               End If
134            End If
135         End Do
136!
137!  The Cs become As and Bs
138!
139         LMTNA = 4
140         Exit
141      End Do
142!
143!  Iteration loop. Each time, the length of the ordered subsets
144!  is doubled.
145!
146      Do
147         If (LMTNA >= NVAL) Exit
148         IWRKF = 0
149         LMTNC = 2 * LMTNC
150!
151!   Loop on merges of A and B into C
152!
153         Do
154            IWRK = IWRKF
155            IWRKD = IWRKF + 1
156            JINDA = IWRKF + LMTNA
157            IWRKF = IWRKF + LMTNC
158            If (IWRKF >= NVAL) Then
159               If (JINDA >= NVAL) Exit
160               IWRKF = NVAL
161            End If
162            IINDA = 1
163            IINDB = JINDA + 1
164!
165!   Shortcut for the case when the max of A is smaller
166!   than the min of B. This line may be activated when the
167!   initial set is already close to sorted.
168!
169!          IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
170!
171!  One steps in the C subset, that we build in the final rank array
172!
173!  Make a copy of the rank array for the merge iteration
174!
175            JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
176!
177            XVALA = XDONT (JWRKT(IINDA))
178            XVALB = XDONT (IRNGT(IINDB))
179!
180            Do
181               IWRK = IWRK + 1
182!
183!  We still have unprocessed values in both A and B
184!
185               If (XVALA > XVALB) Then
186                  IRNGT (IWRK) = IRNGT (IINDB)
187                  IINDB = IINDB + 1
188                  If (IINDB > IWRKF) Then
189!  Only A still with unprocessed values
190                     IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
191                     Exit
192                  End If
193                  XVALB = XDONT (IRNGT(IINDB))
194               Else
195                  IRNGT (IWRK) = JWRKT (IINDA)
196                  IINDA = IINDA + 1
197                  If (IINDA > LMTNA) Exit! Only B still with unprocessed values
198                  XVALA = XDONT (JWRKT(IINDA))
199               End If
200!
201            End Do
202         End Do
203!
204!  The Cs become As and Bs
205!
206         LMTNA = 2 * LMTNA
207      End Do
208!
209      Return
210!
211End Subroutine D_mrgrnk
212
213Subroutine R_mrgrnk (XDONT, IRNGT)
214! __________________________________________________________
215!   MRGRNK = Merge-sort ranking of an array
216!   For performance reasons, the first 2 passes are taken
217!   out of the standard loop, and use dedicated coding.
218! __________________________________________________________
219! _________________________________________________________
220      Real, Dimension (:), Intent (In) :: XDONT
221      Integer*4, Dimension (:), Intent (Out) :: IRNGT
222! __________________________________________________________
223      Real :: XVALA, XVALB
224!
225      Integer, Dimension (SIZE(IRNGT)) :: JWRKT
226      Integer :: LMTNA, LMTNC, IRNG1, IRNG2
227      Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
228!
229      NVAL = Min (SIZE(XDONT), SIZE(IRNGT))
230      Select Case (NVAL)
231      Case (:0)
232         Return
233      Case (1)
234         IRNGT (1) = 1
235         Return
236      Case Default
237         Continue
238      End Select
239!
240!  Fill-in the index array, creating ordered couples
241!
242      Do IIND = 2, NVAL, 2
243         If (XDONT(IIND-1) <= XDONT(IIND)) Then
244            IRNGT (IIND-1) = IIND - 1
245            IRNGT (IIND) = IIND
246         Else
247            IRNGT (IIND-1) = IIND
248            IRNGT (IIND) = IIND - 1
249         End If
250      End Do
251      If (Modulo(NVAL, 2) /= 0) Then
252         IRNGT (NVAL) = NVAL
253      End If
254!
255!  We will now have ordered subsets A - B - A - B - ...
256!  and merge A and B couples into     C   -   C   - ...
257!
258      LMTNA = 2
259      LMTNC = 4
260!
261!  First iteration. The length of the ordered subsets goes from 2 to 4
262!
263      Do
264         If (NVAL <= 2) Exit
265!
266!   Loop on merges of A and B into C
267!
268         Do IWRKD = 0, NVAL - 1, 4
269            If ((IWRKD+4) > NVAL) Then
270               If ((IWRKD+2) >= NVAL) Exit
271!
272!   1 2 3
273!
274               If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
275!
276!   1 3 2
277!
278               If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
279                  IRNG2 = IRNGT (IWRKD+2)
280                  IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
281                  IRNGT (IWRKD+3) = IRNG2
282!
283!   3 1 2
284!
285               Else
286                  IRNG1 = IRNGT (IWRKD+1)
287                  IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
288                  IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
289                  IRNGT (IWRKD+2) = IRNG1
290               End If
291               Exit
292            End If
293!
294!   1 2 3 4
295!
296            If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
297!
298!   1 3 x x
299!
300            If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
301               IRNG2 = IRNGT (IWRKD+2)
302               IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
303               If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
304!   1 3 2 4
305                  IRNGT (IWRKD+3) = IRNG2
306               Else
307!   1 3 4 2
308                  IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
309                  IRNGT (IWRKD+4) = IRNG2
310               End If
311!
312!   3 x x x
313!
314            Else
315               IRNG1 = IRNGT (IWRKD+1)
316               IRNG2 = IRNGT (IWRKD+2)
317               IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
318               If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
319                  IRNGT (IWRKD+2) = IRNG1
320                  If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
321!   3 1 2 4
322                     IRNGT (IWRKD+3) = IRNG2
323                  Else
324!   3 1 4 2
325                     IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
326                     IRNGT (IWRKD+4) = IRNG2
327                  End If
328               Else
329!   3 4 1 2
330                  IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
331                  IRNGT (IWRKD+3) = IRNG1
332                  IRNGT (IWRKD+4) = IRNG2
333               End If
334            End If
335         End Do
336!
337!  The Cs become As and Bs
338!
339         LMTNA = 4
340         Exit
341      End Do
342!
343!  Iteration loop. Each time, the length of the ordered subsets
344!  is doubled.
345!
346      Do
347         If (LMTNA >= NVAL) Exit
348         IWRKF = 0
349         LMTNC = 2 * LMTNC
350!
351!   Loop on merges of A and B into C
352!
353         Do
354            IWRK = IWRKF
355            IWRKD = IWRKF + 1
356            JINDA = IWRKF + LMTNA
357            IWRKF = IWRKF + LMTNC
358            If (IWRKF >= NVAL) Then
359               If (JINDA >= NVAL) Exit
360               IWRKF = NVAL
361            End If
362            IINDA = 1
363            IINDB = JINDA + 1
364!
365!   Shortcut for the case when the max of A is smaller
366!   than the min of B. This line may be activated when the
367!   initial set is already close to sorted.
368!
369!          IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
370!
371!  One steps in the C subset, that we build in the final rank array
372!
373!  Make a copy of the rank array for the merge iteration
374!
375            JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
376!
377            XVALA = XDONT (JWRKT(IINDA))
378            XVALB = XDONT (IRNGT(IINDB))
379!
380            Do
381               IWRK = IWRK + 1
382!
383!  We still have unprocessed values in both A and B
384!
385               If (XVALA > XVALB) Then
386                  IRNGT (IWRK) = IRNGT (IINDB)
387                  IINDB = IINDB + 1
388                  If (IINDB > IWRKF) Then
389!  Only A still with unprocessed values
390                     IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
391                     Exit
392                  End If
393                  XVALB = XDONT (IRNGT(IINDB))
394               Else
395                  IRNGT (IWRK) = JWRKT (IINDA)
396                  IINDA = IINDA + 1
397                  If (IINDA > LMTNA) Exit! Only B still with unprocessed values
398                  XVALA = XDONT (JWRKT(IINDA))
399               End If
400!
401            End Do
402         End Do
403!
404!  The Cs become As and Bs
405!
406         LMTNA = 2 * LMTNA
407      End Do
408!
409      Return
410!
411End Subroutine R_mrgrnk
412Subroutine I_mrgrnk (XDONT, IRNGT)
413! __________________________________________________________
414!   MRGRNK = Merge-sort ranking of an array
415!   For performance reasons, the first 2 passes are taken
416!   out of the standard loop, and use dedicated coding.
417! __________________________________________________________
418! __________________________________________________________
419      Integer, Dimension (:), Intent (In)  :: XDONT
420      Integer*4, Dimension (:), Intent (Out) :: IRNGT
421! __________________________________________________________
422      Integer :: XVALA, XVALB
423!
424      Integer, Dimension (SIZE(IRNGT)) :: JWRKT
425      Integer :: LMTNA, LMTNC, IRNG1, IRNG2
426      Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
427!
428      NVAL = Min (SIZE(XDONT), SIZE(IRNGT))
429      Select Case (NVAL)
430      Case (:0)
431         Return
432      Case (1)
433         IRNGT (1) = 1
434         Return
435      Case Default
436         Continue
437      End Select
438!
439!  Fill-in the index array, creating ordered couples
440!
441      Do IIND = 2, NVAL, 2
442         If (XDONT(IIND-1) <= XDONT(IIND)) Then
443            IRNGT (IIND-1) = IIND - 1
444            IRNGT (IIND) = IIND
445         Else
446            IRNGT (IIND-1) = IIND
447            IRNGT (IIND) = IIND - 1
448         End If
449      End Do
450      If (Modulo(NVAL, 2) /= 0) Then
451         IRNGT (NVAL) = NVAL
452      End If
453!
454!  We will now have ordered subsets A - B - A - B - ...
455!  and merge A and B couples into     C   -   C   - ...
456!
457      LMTNA = 2
458      LMTNC = 4
459!
460!  First iteration. The length of the ordered subsets goes from 2 to 4
461!
462      Do
463         If (NVAL <= 2) Exit
464!
465!   Loop on merges of A and B into C
466!
467         Do IWRKD = 0, NVAL - 1, 4
468            If ((IWRKD+4) > NVAL) Then
469               If ((IWRKD+2) >= NVAL) Exit
470!
471!   1 2 3
472!
473               If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
474!
475!   1 3 2
476!
477               If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
478                  IRNG2 = IRNGT (IWRKD+2)
479                  IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
480                  IRNGT (IWRKD+3) = IRNG2
481!
482!   3 1 2
483!
484               Else
485                  IRNG1 = IRNGT (IWRKD+1)
486                  IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
487                  IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
488                  IRNGT (IWRKD+2) = IRNG1
489               End If
490               Exit
491            End If
492!
493!   1 2 3 4
494!
495            If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
496!
497!   1 3 x x
498!
499            If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
500               IRNG2 = IRNGT (IWRKD+2)
501               IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
502               If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
503!   1 3 2 4
504                  IRNGT (IWRKD+3) = IRNG2
505               Else
506!   1 3 4 2
507                  IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
508                  IRNGT (IWRKD+4) = IRNG2
509               End If
510!
511!   3 x x x
512!
513            Else
514               IRNG1 = IRNGT (IWRKD+1)
515               IRNG2 = IRNGT (IWRKD+2)
516               IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
517               If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
518                  IRNGT (IWRKD+2) = IRNG1
519                  If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
520!   3 1 2 4
521                     IRNGT (IWRKD+3) = IRNG2
522                  Else
523!   3 1 4 2
524                     IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
525                     IRNGT (IWRKD+4) = IRNG2
526                  End If
527               Else
528!   3 4 1 2
529                  IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
530                  IRNGT (IWRKD+3) = IRNG1
531                  IRNGT (IWRKD+4) = IRNG2
532               End If
533            End If
534         End Do
535!
536!  The Cs become As and Bs
537!
538         LMTNA = 4
539         Exit
540      End Do
541!
542!  Iteration loop. Each time, the length of the ordered subsets
543!  is doubled.
544!
545      Do
546         If (LMTNA >= NVAL) Exit
547         IWRKF = 0
548         LMTNC = 2 * LMTNC
549!
550!   Loop on merges of A and B into C
551!
552         Do
553            IWRK = IWRKF
554            IWRKD = IWRKF + 1
555            JINDA = IWRKF + LMTNA
556            IWRKF = IWRKF + LMTNC
557            If (IWRKF >= NVAL) Then
558               If (JINDA >= NVAL) Exit
559               IWRKF = NVAL
560            End If
561            IINDA = 1
562            IINDB = JINDA + 1
563!
564!   Shortcut for the case when the max of A is smaller
565!   than the min of B. This line may be activated when the
566!   initial set is already close to sorted.
567!
568!          IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
569!
570!  One steps in the C subset, that we build in the final rank array
571!
572!  Make a copy of the rank array for the merge iteration
573!
574            JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
575!
576            XVALA = XDONT (JWRKT(IINDA))
577            XVALB = XDONT (IRNGT(IINDB))
578!
579            Do
580               IWRK = IWRK + 1
581!
582!  We still have unprocessed values in both A and B
583!
584               If (XVALA > XVALB) Then
585                  IRNGT (IWRK) = IRNGT (IINDB)
586                  IINDB = IINDB + 1
587                  If (IINDB > IWRKF) Then
588!  Only A still with unprocessed values
589                     IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
590                     Exit
591                  End If
592                  XVALB = XDONT (IRNGT(IINDB))
593               Else
594                  IRNGT (IWRK) = JWRKT (IINDA)
595                  IINDA = IINDA + 1
596                  If (IINDA > LMTNA) Exit! Only B still with unprocessed values
597                  XVALA = XDONT (JWRKT(IINDA))
598               End If
599!
600            End Do
601         End Do
602!
603!  The Cs become As and Bs
604!
605         LMTNA = 2 * LMTNA
606      End Do
607!
608      Return
609!
610End Subroutine I_mrgrnk
611end module m_mrgrnk
Note: See TracBrowser for help on using the repository browser.