source: LMDZ6/trunk/libf/phylmd/cospv2/mrgrnk.F90 @ 3825

Last change on this file since 3825 was 3491, checked in by idelkadi, 6 years ago

Integration of version 2 of the COSP simulator in LMDZ
This line, and those below, will be ignored--

M makegcm
M makelmdz
M makelmdz_fcm
M libf/phylmd/physiq_mod.F90
A libf/phylmd/cospv2
A libf/phylmd/cospv2/mo_rng.F90
A libf/phylmd/cospv2/quickbeam_optics.F90
A libf/phylmd/cospv2/cosp_cloudsat_interface.F90
A libf/phylmd/cospv2/cosp_config.F90
A libf/phylmd/cospv2/lidar_simulator.F90
A libf/phylmd/cospv2/prec_scops.F90
A libf/phylmd/cospv2/mrgrnk.F90
A libf/phylmd/cospv2/lmdz_cosp_read_outputkeys.F90
A libf/phylmd/cospv2/cosp_atlid_interface.F90
A libf/phylmd/cospv2/lmdz_cosp_subsample_and_optics_mod.F90
A libf/phylmd/cospv2/cosp_math_constants.F90
A libf/phylmd/cospv2/MISR_simulator.F90
A libf/phylmd/cospv2/modis_simulator.F90
A libf/phylmd/cospv2/math_lib.F90
A libf/phylmd/cospv2/cosp_grLidar532_interface.F90
A libf/phylmd/cospv2/cosp_errorHandling.F90
A libf/phylmd/cospv2/cosp_stats.F90
A libf/phylmd/cospv2/lmdz_cosp_output_write_mod.F90
A libf/phylmd/cospv2/cosp_utils.F90
A libf/phylmd/cospv2/cosp_optics.F90
A libf/phylmd/cospv2/icarus.F90
A libf/phylmd/cospv2/scops.F90
A libf/phylmd/cospv2/optics_lib.F90
A libf/phylmd/cospv2/cosp_kinds.F90
A libf/phylmd/cospv2/cosp_calipso_interface.F90
A libf/phylmd/cospv2/quickbeam.F90
A libf/phylmd/cospv2/parasol.F90
A libf/phylmd/cospv2/cosp_phys_constants.F90
A libf/phylmd/cospv2/cosp.F90
A libf/phylmd/cospv2/array_lib.F90
A libf/phylmd/cospv2/cosp_isccp_interface.F90
A libf/phylmd/cospv2/cosp_parasol_interface.F90
A libf/phylmd/cospv2/lmdz_cosp_construct_destroy_mod.F90
A libf/phylmd/cospv2/lmdz_cosp_output_mod.F90
A libf/phylmd/cospv2/lmdz_cosp_interface.F90
A libf/phylmd/cospv2/cosp_misr_interface.F90
A libf/phylmd/cospv2/cosp_modis_interface.F90
File size: 19.7 KB
Line 
1! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2! Copyright (c) 2015, Regents of the University of Colorado
3! All rights reserved.
4!
5! Redistribution and use in source and binary forms, with or without modification, are
6! permitted provided that the following conditions are met:
7!
8! 1. Redistributions of source code must retain the above copyright notice, this list of
9!    conditions and the following disclaimer.
10!
11! 2. Redistributions in binary form must reproduce the above copyright notice, this list
12!    of conditions and the following disclaimer in the documentation and/or other
13!    materials provided with the distribution.
14!
15! 3. Neither the name of the copyright holder nor the names of its contributors may be
16!    used to endorse or promote products derived from this software without specific prior
17!    written permission.
18!
19! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
20! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
21! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
22! THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
23! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
24! OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28!
29! History:
30! May 2015:  Dustin Swales    - Modified for COSPv2.0
31!
32! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
33Module m_mrgrnk
34  USE COSP_KINDS,          ONLY: wp
35  Integer, Parameter :: kdp = selected_real_kind(15)
36  public :: mrgrnk
37  private :: kdp
38  private :: R_mrgrnk, I_mrgrnk, D_mrgrnk
39
40  interface mrgrnk
41!     module procedure D_mrgrnk, R_mrgrnk, I_mrgrnk
42     module procedure R_mrgrnk, I_mrgrnk
43     
44  end interface
45contains
46 
47  Subroutine D_mrgrnk (XDONT, IRNGT)
48    ! __________________________________________________________
49    ! MRGRNK = Merge-sort ranking of an array
50    ! For performance reasons, the first 2 passes are taken
51    ! out of the standard loop, and use dedicated coding.
52    ! __________________________________________________________
53    ! __________________________________________________________
54    Real (wp), Dimension (:), Intent (In) :: XDONT
55    Integer, Dimension (:), Intent (Out) :: IRNGT
56    ! __________________________________________________________
57    Real (wp) :: XVALA, XVALB
58    !
59    Integer, Dimension (SIZE(IRNGT)) :: JWRKT
60    Integer :: LMTNA, LMTNC, IRNG1, IRNG2
61    Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
62    !
63    NVAL = Min (SIZE(XDONT), SIZE(IRNGT))
64    Select Case (NVAL)
65    Case (:0)
66       Return
67    Case (1)
68       IRNGT (1) = 1
69       Return
70    Case Default
71       Continue
72    End Select
73    !
74    ! Fill-in the index array, creating ordered couples
75    !
76    Do IIND = 2, NVAL, 2
77       If (XDONT(IIND-1) <= XDONT(IIND)) Then
78          IRNGT (IIND-1) = IIND - 1
79          IRNGT (IIND) = IIND
80       Else
81          IRNGT (IIND-1) = IIND
82          IRNGT (IIND) = IIND - 1
83       End If
84    End Do
85    If (Modulo(NVAL, 2) /= 0) Then
86       IRNGT (NVAL) = NVAL
87    End If
88    !
89    ! We will now have ordered subsets A - B - A - B - ...
90    ! and merge A and B couples into C - C - ...
91    !
92    LMTNA = 2
93    LMTNC = 4
94    !
95    ! First iteration. The length of the ordered subsets goes from 2 to 4
96    !
97    Do
98       If (NVAL <= 2) Exit
99       !
100       ! Loop on merges of A and B into C
101       !
102       Do IWRKD = 0, NVAL - 1, 4
103          If ((IWRKD+4) > NVAL) Then
104             If ((IWRKD+2) >= NVAL) Exit
105             !
106             ! 1 2 3
107             !
108             If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
109             !
110             ! 1 3 2
111             !
112             If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
113                IRNG2 = IRNGT (IWRKD+2)
114                IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
115                IRNGT (IWRKD+3) = IRNG2
116                !
117                ! 3 1 2
118                !
119             Else
120                IRNG1 = IRNGT (IWRKD+1)
121                IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
122                IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
123                IRNGT (IWRKD+2) = IRNG1
124             End If
125             Exit
126          End If
127          !
128          ! 1 2 3 4
129          !
130          If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
131          !
132          ! 1 3 x x
133          !
134          If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
135             IRNG2 = IRNGT (IWRKD+2)
136             IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
137             If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
138                ! 1 3 2 4
139                IRNGT (IWRKD+3) = IRNG2
140             Else
141                ! 1 3 4 2
142                IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
143                IRNGT (IWRKD+4) = IRNG2
144             End If
145             !
146             ! 3 x x x
147             !
148          Else
149             IRNG1 = IRNGT (IWRKD+1)
150             IRNG2 = IRNGT (IWRKD+2)
151             IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
152             If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
153                IRNGT (IWRKD+2) = IRNG1
154                If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
155                   ! 3 1 2 4
156                   IRNGT (IWRKD+3) = IRNG2
157                Else
158                   ! 3 1 4 2
159                   IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
160                   IRNGT (IWRKD+4) = IRNG2
161                End If
162             Else
163                ! 3 4 1 2
164                IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
165                IRNGT (IWRKD+3) = IRNG1
166                IRNGT (IWRKD+4) = IRNG2
167             End If
168          End If
169       End Do
170       !
171       ! The Cs become As and Bs
172       !
173       LMTNA = 4
174       Exit
175    End Do
176    !
177    ! Iteration loop. Each time, the length of the ordered subsets
178    ! is doubled.
179    !
180    Do
181       If (LMTNA >= NVAL) Exit
182       IWRKF = 0
183       LMTNC = 2 * LMTNC
184       !
185       ! Loop on merges of A and B into C
186       !
187       Do
188          IWRK = IWRKF
189          IWRKD = IWRKF + 1
190          JINDA = IWRKF + LMTNA
191          IWRKF = IWRKF + LMTNC
192          If (IWRKF >= NVAL) Then
193             If (JINDA >= NVAL) Exit
194             IWRKF = NVAL
195          End If
196          IINDA = 1
197          IINDB = JINDA + 1
198          !
199          ! Shortcut for the case when the max of A is smaller
200          ! than the min of B. This line may be activated when the
201          ! initial set is already close to sorted.
202          !
203          ! IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
204          !
205          ! One steps in the C subset, that we build in the final rank array
206          !
207          ! Make a copy of the rank array for the merge iteration
208          !
209          JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
210          !
211          XVALA = XDONT (JWRKT(IINDA))
212          XVALB = XDONT (IRNGT(IINDB))
213          !
214          Do
215             IWRK = IWRK + 1
216             !
217             ! We still have unprocessed values in both A and B
218             !
219             If (XVALA > XVALB) Then
220                IRNGT (IWRK) = IRNGT (IINDB)
221                IINDB = IINDB + 1
222                If (IINDB > IWRKF) Then
223                   ! Only A still with unprocessed values
224                   IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
225                   Exit
226                End If
227                XVALB = XDONT (IRNGT(IINDB))
228             Else
229                IRNGT (IWRK) = JWRKT (IINDA)
230                IINDA = IINDA + 1
231                If (IINDA > LMTNA) Exit! Only B still with unprocessed values
232                XVALA = XDONT (JWRKT(IINDA))
233             End If
234             !
235          End Do
236       End Do
237       !
238       ! The Cs become As and Bs
239       !
240       LMTNA = 2 * LMTNA
241    End Do
242    !
243    Return
244    !
245  End Subroutine D_mrgrnk
246 
247  Subroutine R_mrgrnk (XDONT, IRNGT)
248    ! __________________________________________________________
249    ! MRGRNK = Merge-sort ranking of an array
250    ! For performance reasons, the first 2 passes are taken
251    ! out of the standard loop, and use dedicated coding.
252    ! __________________________________________________________
253    ! _________________________________________________________
254    Real(wp), Dimension (:), Intent (In) :: XDONT
255    Integer, Dimension (:), Intent (Out) :: IRNGT
256    ! __________________________________________________________
257    Real(wp) :: XVALA, XVALB
258    !
259    Integer, Dimension (SIZE(IRNGT)) :: JWRKT
260    Integer :: LMTNA, LMTNC, IRNG1, IRNG2
261    Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
262    !
263    NVAL = Min (SIZE(XDONT), SIZE(IRNGT))
264    Select Case (NVAL)
265    Case (:0)
266       Return
267    Case (1)
268       IRNGT (1) = 1
269       Return
270    Case Default
271       Continue
272    End Select
273    !
274    ! Fill-in the index array, creating ordered couples
275    !
276    Do IIND = 2, NVAL, 2
277       If (XDONT(IIND-1) <= XDONT(IIND)) Then
278          IRNGT (IIND-1) = IIND - 1
279          IRNGT (IIND) = IIND
280       Else
281          IRNGT (IIND-1) = IIND
282          IRNGT (IIND) = IIND - 1
283       End If
284    End Do
285    If (Modulo(NVAL, 2) /= 0) Then
286       IRNGT (NVAL) = NVAL
287    End If
288    !
289    ! We will now have ordered subsets A - B - A - B - ...
290    ! and merge A and B couples into C - C - ...
291    !
292    LMTNA = 2
293    LMTNC = 4
294    !
295    ! First iteration. The length of the ordered subsets goes from 2 to 4
296    !
297    Do
298       If (NVAL <= 2) Exit
299       !
300       ! Loop on merges of A and B into C
301       !
302       Do IWRKD = 0, NVAL - 1, 4
303          If ((IWRKD+4) > NVAL) Then
304             If ((IWRKD+2) >= NVAL) Exit
305             !
306             ! 1 2 3
307             !
308             If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
309             !
310             ! 1 3 2
311             !
312             If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
313                IRNG2 = IRNGT (IWRKD+2)
314                IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
315                IRNGT (IWRKD+3) = IRNG2
316                !
317                ! 3 1 2
318                !
319             Else
320                IRNG1 = IRNGT (IWRKD+1)
321                IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
322                IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
323                IRNGT (IWRKD+2) = IRNG1
324             End If
325             Exit
326          End If
327          !
328          ! 1 2 3 4
329          !
330          If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
331          !
332          ! 1 3 x x
333          !
334          If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
335             IRNG2 = IRNGT (IWRKD+2)
336             IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
337             If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
338                ! 1 3 2 4
339                IRNGT (IWRKD+3) = IRNG2
340             Else
341                ! 1 3 4 2
342                IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
343                IRNGT (IWRKD+4) = IRNG2
344             End If
345             !
346             ! 3 x x x
347             !
348          Else
349             IRNG1 = IRNGT (IWRKD+1)
350             IRNG2 = IRNGT (IWRKD+2)
351             IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
352             If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
353                IRNGT (IWRKD+2) = IRNG1
354                If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
355                   ! 3 1 2 4
356                   IRNGT (IWRKD+3) = IRNG2
357                Else
358                   ! 3 1 4 2
359                   IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
360                   IRNGT (IWRKD+4) = IRNG2
361                End If
362             Else
363                ! 3 4 1 2
364                IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
365                IRNGT (IWRKD+3) = IRNG1
366                IRNGT (IWRKD+4) = IRNG2
367             End If
368          End If
369       End Do
370       !
371       ! The Cs become As and Bs
372       !
373       LMTNA = 4
374       Exit
375    End Do
376    !
377    ! Iteration loop. Each time, the length of the ordered subsets
378    ! is doubled.
379    !
380    Do
381       If (LMTNA >= NVAL) Exit
382       IWRKF = 0
383       LMTNC = 2 * LMTNC
384       !
385       ! Loop on merges of A and B into C
386       !
387       Do
388          IWRK = IWRKF
389          IWRKD = IWRKF + 1
390          JINDA = IWRKF + LMTNA
391          IWRKF = IWRKF + LMTNC
392          If (IWRKF >= NVAL) Then
393             If (JINDA >= NVAL) Exit
394             IWRKF = NVAL
395          End If
396          IINDA = 1
397          IINDB = JINDA + 1
398          !
399          ! Shortcut for the case when the max of A is smaller
400          ! than the min of B. This line may be activated when the
401          ! initial set is already close to sorted.
402          !
403          ! IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
404          !
405          ! One steps in the C subset, that we build in the final rank array
406          !
407          ! Make a copy of the rank array for the merge iteration
408          !
409          JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
410          !
411          XVALA = XDONT (JWRKT(IINDA))
412          XVALB = XDONT (IRNGT(IINDB))
413          !
414          Do
415             IWRK = IWRK + 1
416             !
417             ! We still have unprocessed values in both A and B
418             !
419             If (XVALA > XVALB) Then
420                IRNGT (IWRK) = IRNGT (IINDB)
421                IINDB = IINDB + 1
422                If (IINDB > IWRKF) Then
423                   ! Only A still with unprocessed values
424                   IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
425                   Exit
426                End If
427                XVALB = XDONT (IRNGT(IINDB))
428             Else
429                IRNGT (IWRK) = JWRKT (IINDA)
430                IINDA = IINDA + 1
431                If (IINDA > LMTNA) Exit! Only B still with unprocessed values
432                XVALA = XDONT (JWRKT(IINDA))
433             End If
434             !
435          End Do
436       End Do
437       !
438       ! The Cs become As and Bs
439       !
440       LMTNA = 2 * LMTNA
441    End Do
442    !
443    Return
444    !
445  End Subroutine R_mrgrnk
446  Subroutine I_mrgrnk (XDONT, IRNGT)
447    ! __________________________________________________________
448    ! MRGRNK = Merge-sort ranking of an array
449    ! For performance reasons, the first 2 passes are taken
450    ! out of the standard loop, and use dedicated coding.
451    ! __________________________________________________________
452    ! __________________________________________________________
453    Integer, Dimension (:), Intent (In) :: XDONT
454    Integer, Dimension (:), Intent (Out) :: IRNGT
455    ! __________________________________________________________
456    Integer :: XVALA, XVALB
457    !
458    Integer, Dimension (SIZE(IRNGT)) :: JWRKT
459    Integer :: LMTNA, LMTNC, IRNG1, IRNG2
460    Integer :: NVAL, IIND, IWRKD, IWRK, IWRKF, JINDA, IINDA, IINDB
461    !
462    NVAL = Min (SIZE(XDONT), SIZE(IRNGT))
463    Select Case (NVAL)
464    Case (:0)
465       Return
466    Case (1)
467       IRNGT (1) = 1
468       Return
469    Case Default
470       Continue
471    End Select
472    !
473    ! Fill-in the index array, creating ordered couples
474    !
475    Do IIND = 2, NVAL, 2
476       If (XDONT(IIND-1) <= XDONT(IIND)) Then
477          IRNGT (IIND-1) = IIND - 1
478          IRNGT (IIND) = IIND
479       Else
480          IRNGT (IIND-1) = IIND
481          IRNGT (IIND) = IIND - 1
482       End If
483    End Do
484    If (Modulo(NVAL, 2) /= 0) Then
485       IRNGT (NVAL) = NVAL
486    End If
487    !
488    ! We will now have ordered subsets A - B - A - B - ...
489    ! and merge A and B couples into C - C - ...
490    !
491    LMTNA = 2
492    LMTNC = 4
493    !
494    ! First iteration. The length of the ordered subsets goes from 2 to 4
495    !
496    Do
497       If (NVAL <= 2) Exit
498       !
499       ! Loop on merges of A and B into C
500       !
501       Do IWRKD = 0, NVAL - 1, 4
502          If ((IWRKD+4) > NVAL) Then
503             If ((IWRKD+2) >= NVAL) Exit
504             !
505             ! 1 2 3
506             !
507             If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Exit
508             !
509             ! 1 3 2
510             !
511             If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
512                IRNG2 = IRNGT (IWRKD+2)
513                IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
514                IRNGT (IWRKD+3) = IRNG2
515                !
516                ! 3 1 2
517                !
518             Else
519                IRNG1 = IRNGT (IWRKD+1)
520                IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
521                IRNGT (IWRKD+3) = IRNGT (IWRKD+2)
522                IRNGT (IWRKD+2) = IRNG1
523             End If
524             Exit
525          End If
526          !
527          ! 1 2 3 4
528          !
529          If (XDONT(IRNGT(IWRKD+2)) <= XDONT(IRNGT(IWRKD+3))) Cycle
530          !
531          ! 1 3 x x
532          !
533          If (XDONT(IRNGT(IWRKD+1)) <= XDONT(IRNGT(IWRKD+3))) Then
534             IRNG2 = IRNGT (IWRKD+2)
535             IRNGT (IWRKD+2) = IRNGT (IWRKD+3)
536             If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
537                ! 1 3 2 4
538                IRNGT (IWRKD+3) = IRNG2
539             Else
540                ! 1 3 4 2
541                IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
542                IRNGT (IWRKD+4) = IRNG2
543             End If
544             !
545             ! 3 x x x
546             !
547          Else
548             IRNG1 = IRNGT (IWRKD+1)
549             IRNG2 = IRNGT (IWRKD+2)
550             IRNGT (IWRKD+1) = IRNGT (IWRKD+3)
551             If (XDONT(IRNG1) <= XDONT(IRNGT(IWRKD+4))) Then
552                IRNGT (IWRKD+2) = IRNG1
553                If (XDONT(IRNG2) <= XDONT(IRNGT(IWRKD+4))) Then
554                   ! 3 1 2 4
555                   IRNGT (IWRKD+3) = IRNG2
556                Else
557                   ! 3 1 4 2
558                   IRNGT (IWRKD+3) = IRNGT (IWRKD+4)
559                   IRNGT (IWRKD+4) = IRNG2
560                End If
561             Else
562                ! 3 4 1 2
563                IRNGT (IWRKD+2) = IRNGT (IWRKD+4)
564                IRNGT (IWRKD+3) = IRNG1
565                IRNGT (IWRKD+4) = IRNG2
566             End If
567          End If
568       End Do
569       !
570       ! The Cs become As and Bs
571       !
572       LMTNA = 4
573       Exit
574    End Do
575    !
576    ! Iteration loop. Each time, the length of the ordered subsets
577    ! is doubled.
578    !
579    Do
580       If (LMTNA >= NVAL) Exit
581       IWRKF = 0
582       LMTNC = 2 * LMTNC
583       !
584       ! Loop on merges of A and B into C
585       !
586       Do
587          IWRK = IWRKF
588          IWRKD = IWRKF + 1
589          JINDA = IWRKF + LMTNA
590          IWRKF = IWRKF + LMTNC
591          If (IWRKF >= NVAL) Then
592             If (JINDA >= NVAL) Exit
593             IWRKF = NVAL
594          End If
595          IINDA = 1
596          IINDB = JINDA + 1
597          !
598          ! Shortcut for the case when the max of A is smaller
599          ! than the min of B. This line may be activated when the
600          ! initial set is already close to sorted.
601          !
602          ! IF (XDONT(IRNGT(JINDA)) <= XDONT(IRNGT(IINDB))) CYCLE
603          !
604          ! One steps in the C subset, that we build in the final rank array
605          !
606          ! Make a copy of the rank array for the merge iteration
607          !
608          JWRKT (1:LMTNA) = IRNGT (IWRKD:JINDA)
609          !
610          XVALA = XDONT (JWRKT(IINDA))
611          XVALB = XDONT (IRNGT(IINDB))
612          !
613          Do
614             IWRK = IWRK + 1
615             !
616             ! We still have unprocessed values in both A and B
617             !
618             If (XVALA > XVALB) Then
619                IRNGT (IWRK) = IRNGT (IINDB)
620                IINDB = IINDB + 1
621                If (IINDB > IWRKF) Then
622                   ! Only A still with unprocessed values
623                   IRNGT (IWRK+1:IWRKF) = JWRKT (IINDA:LMTNA)
624                   Exit
625                End If
626                XVALB = XDONT (IRNGT(IINDB))
627             Else
628                IRNGT (IWRK) = JWRKT (IINDA)
629                IINDA = IINDA + 1
630                If (IINDA > LMTNA) Exit! Only B still with unprocessed values
631                XVALA = XDONT (JWRKT(IINDA))
632             End If
633             !
634          End Do
635       End Do
636       !
637       ! The Cs become As and Bs
638       !
639       LMTNA = 2 * LMTNA
640    End Do
641    !
642    Return
643    !
644  End Subroutine I_mrgrnk
645end module m_mrgrnk
Note: See TracBrowser for help on using the repository browser.