Ignore:
Timestamp:
Jun 11, 2020 8:51:48 AM (11 months ago)
Author:
raasch
Message:

files re-formatted to follow the PALM coding standard

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/temperton_fft_mod.f90

    r4370 r4559  
    11!> @file temperton_fft_mod.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
     3! This file is part of the PALM model system.
     4!
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
     8!
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
     15!
     16! Copyright 1997-2020 Leibniz Universitaet Hannover
     17!--------------------------------------------------------------------------------------------------!
     18!
    319!
    420! Current revisions:
     
    925! -----------------
    1026! $Id$
     27! File re-formatted to follow the PALM coding standard
     28!
     29! 4370 2020-01-10 14:00:44Z raasch
    1130! unused variables removed
    12 ! 
     31!
    1332! 4366 2020-01-09 08:12:43Z raasch
    1433! vectorized routines added
    15 ! 
     34!
    1635! 4182 2019-08-22 15:20:23Z scharf
    1736! Corrected "Former revisions" section
    18 ! 
     37!
    1938! 3725 2019-02-07 10:11:02Z raasch
    2039! GOTO statements replaced, file re-formatted corresponding to coding standards
     
    2746! ------------
    2847!> Fast Fourier transformation developed by Clive Temperton, ECMWF.
    29 !------------------------------------------------------------------------------!
     48!--------------------------------------------------------------------------------------------------!
    3049 MODULE temperton_fft
    3150
     
    5776CONTAINS
    5877
    59 !------------------------------------------------------------------------------!
     78!--------------------------------------------------------------------------------------------------!
    6079! Description:
    6180! ------------
    6281!> Calls Fortran-versions of fft's.
    63 !> 
     82!>
    6483!> Method:
    65 !> 
    66 !> Subroutine 'fft991cy' - multiple fast real periodic transform
    67 !> supersedes previous routine 'fft991cy'.
    68 !> 
    69 !> Real transform of length n performed by removing redundant
    70 !> operations from complex transform of length n.
    71 !> 
     84!>
     85!> Subroutine 'fft991cy' - multiple fast real periodic transform supersedes previous routine
     86!> 'fft991cy'.
     87!>
     88!> Real transform of length n performed by removing redundant operations from complex transform of
     89!> length n.
     90!>
    7291!> a       is the array containing input & output data.
    7392!> work    is an area of size (n+1)*min(lot,nfft).
     
    81100!> isign = +1 for transform from spectral to gridpoint
    82101!>       = -1 for transform from gridpoint to spectral.
    83 !> 
     102!>
    84103!> ordering of coefficients:
    85104!> a(0),b(0),a(1),b(1),a(2),b(2),.,a(n/2),b(n/2)
    86105!> where b(0)=b(n/2)=0; (n+2) locations required.
    87 !> 
     106!>
    88107!> ordering of data:
    89108!> x(0),x(1),x(2),.,x(n-1), 0 , 0 ; (n+2) locations required.
    90 !>
    91 !> Vectorization is achieved on cray by doing the transforms
    92 !> in parallel.
    93 !>
     109!>
     110!> Vectorization is achieved on cray by doing the transforms in parallel.
     111!>
    94112!> n must be composed of factors 2,3 & 5 but does not have to be even.
    95 !> 
     113!>
    96114!> definition of transforms:
    97 !> 
     115!>
    98116!> isign=+1: x(j)=sum(k=0,.,n-1)(c(k)*exp(2*i*j*k*pi/n))
    99117!> where c(k)=a(k)+i*b(k) and c(n-k)=a(k)-i*b(k)
    100 !> 
     118!>
    101119!> isign=-1: a(k)=(1/n)*sum(j=0,.,n-1)(x(j)*cos(2*j*k*pi/n))
    102120!> b(k)=-(1/n)*sum(j=0,.,n-1)(x(j)*sin(2*j*k*pi/n))
    103 !> 
     121!>
    104122!> calls fortran-versions of fft's  !!!
    105123!> dimension a(n),work(n),trigs(n),ifax(1)
    106 !------------------------------------------------------------------------------!
     124!--------------------------------------------------------------------------------------------------!
    107125 SUBROUTINE fft991cy( a, work, trigs, ifax, inc, jump, n, lot, isign )
    108126
     
    111129    IMPLICIT NONE
    112130
    113     INTEGER(iwp) ::  inc   !<
    114     INTEGER(iwp) ::  isign !<
    115     INTEGER(iwp) ::  jump  !<
    116     INTEGER(iwp) ::  lot   !<
    117     INTEGER(iwp) ::  n     !<
    118 
    119     REAL(wp)     ::  a(*)     !<
    120     REAL(wp)     ::  trigs(*) !<
    121     REAL(wp)     ::  work(*)  !<
    122     INTEGER(iwp) ::  ifax(*)  !<
    123 
    124     INTEGER(iwp) ::  i      !<
    125     INTEGER(iwp) ::  ia     !<
    126     INTEGER(iwp) ::  ibase  !<
    127     INTEGER(iwp) ::  ierr   !<
    128     INTEGER(iwp) ::  ifac   !<
    129     INTEGER(iwp) ::  igo    !<
    130     INTEGER(iwp) ::  ii     !<
    131     INTEGER(iwp) ::  istart !<
    132     INTEGER(iwp) ::  ix     !<
    133     INTEGER(iwp) ::  iz     !<
    134     INTEGER(iwp) ::  j      !<
    135     INTEGER(iwp) ::  jbase  !<
    136     INTEGER(iwp) ::  jj     !<
    137     INTEGER(iwp) ::  k      !<
    138     INTEGER(iwp) ::  la     !<
    139     INTEGER(iwp) ::  nb     !<
    140     INTEGER(iwp) ::  nblox  !<
    141     INTEGER(iwp) ::  nfax   !<
    142     INTEGER(iwp) ::  nvex   !<
    143     INTEGER(iwp) ::  nx     !<
     131    INTEGER(iwp) ::  inc    !<
     132    INTEGER(iwp) ::  isign  !<
     133    INTEGER(iwp) ::  jump   !<
     134    INTEGER(iwp) ::  lot    !<
     135    INTEGER(iwp) ::  n      !<
     136
     137    INTEGER(iwp) ::  i        !<
     138    INTEGER(iwp) ::  ia       !<
     139    INTEGER(iwp) ::  ibase    !<
     140    INTEGER(iwp) ::  ierr     !<
     141    INTEGER(iwp) ::  ifac     !<
     142    INTEGER(iwp) ::  igo      !<
     143    INTEGER(iwp) ::  ii       !<
     144    INTEGER(iwp) ::  ifax(*)  !<
     145    INTEGER(iwp) ::  istart   !<
     146    INTEGER(iwp) ::  ix       !<
     147    INTEGER(iwp) ::  iz       !<
     148    INTEGER(iwp) ::  j        !<
     149    INTEGER(iwp) ::  jbase    !<
     150    INTEGER(iwp) ::  jj       !<
     151    INTEGER(iwp) ::  k        !<
     152    INTEGER(iwp) ::  la       !<
     153    INTEGER(iwp) ::  nb       !<
     154    INTEGER(iwp) ::  nblox    !<
     155    INTEGER(iwp) ::  nfax     !<
     156    INTEGER(iwp) ::  nvex     !<
     157    INTEGER(iwp) ::  nx       !<
     158
     159    REAL(wp)     ::  a(*)      !<
     160    REAL(wp)     ::  trigs(*)  !<
     161    REAL(wp)     ::  work(*)   !<
    144162
    145163
     
    147165    nfax = ifax(1)
    148166    nx   = n + 1
    149     IF ( MOD(n,2) == 1 )  nx = n
     167    IF ( MOD( n, 2) == 1 )  nx = n
    150168    nblox = 1 + (lot-1) / nfft
    151     nvex = lot - (nblox-1) * nfft
     169    nvex = lot - ( nblox - 1 ) * nfft
    152170
    153171
     
    160178
    161179          ia = istart
    162           i = istart
     180          i  = istart
    163181          !DIR$ IVDEP
    164182          !CDIR NODEP
     
    166184          DO  j = 1, nvex
    167185             a(i+inc) = 0.5_wp * a(i)
    168              i = i + jump
     186             i        = i + jump
    169187          ENDDO
    170188          IF ( MOD(n,2) /= 1 )  THEN
    171189             i = istart + n * inc
    172              DO  j = 1, nvex
     190             DO  j   = 1, nvex
    173191                a(i) = 0.5_wp * a(i)
    174                 i = i + jump
     192                i    = i + jump
    175193             ENDDO
    176194          ENDIF
     
    211229!
    212230!--       If necessary, copy results back to a
    213           IF ( MOD(nfax,2) /= 0 )  THEN
     231          IF ( MOD( nfax, 2 ) /= 0 )  THEN
    214232             ibase = 1
    215233             jbase = ia
     
    219237                DO  ii = 1, n
    220238                   a(j) = work(i)
    221                    i = i + 1
    222                    j = j + inc
     239                   i    = i + 1
     240                   j    = j + inc
    223241                ENDDO
    224242                ibase = ibase + nx
     
    228246!
    229247!--       Fill in zeros at end
    230           ix = istart + n*inc
     248          ix = istart + n * inc
    231249          !DIR$ IVDEP
    232250          !CDIR NODEP
    233251          !OCL NOVREC
    234252          DO  j = 1, nvex
    235              a(ix) = 0.0_wp
     253             a(ix)     = 0.0_wp
    236254             a(ix+inc) = 0.0_wp
    237              ix = ix + jump
     255             ix        = ix + jump
    238256          ENDDO
    239           istart = istart + nvex*jump
     257          istart = istart + nvex * jump
    240258          nvex = nfft
    241259
     
    254272
    255273             ifac = ifax(nfax+2-k)
    256              la = la / ifac
     274             la   = la / ifac
    257275             ierr = -1
    258276             IF ( igo /= -1 )  THEN
     
    291309                DO  ii = 1, n
    292310                   a(j) = work(i)
    293                    i = i + 1
    294                    j = j + inc
     311                   i    = i + 1
     312                   j    = j + inc
    295313                ENDDO
    296314                ibase = ibase + nx
     
    305323          !OCL NOVREC
    306324          DO  j = 1, nvex
    307              a(ix) = a(ix+inc)
     325             a(ix)     = a(ix+inc)
    308326             a(ix+inc) = 0.0_wp
    309              ix = ix + jump
     327             ix        = ix + jump
    310328          ENDDO
    311329          IF ( MOD(n,2) /= 1 )  THEN
     
    313331             DO  j = 1, nvex
    314332                a(iz) = 0.0_wp
    315                 iz = iz + jump
     333                iz    = iz + jump
    316334             ENDDO
    317335          ENDIF
     
    325343 END SUBROUTINE fft991cy
    326344
    327 !------------------------------------------------------------------------------!
     345!--------------------------------------------------------------------------------------------------!
    328346! Description:
    329347! ------------
    330 !> Performs one pass through data as part of
    331 !> multiple real fft (fourier analysis) routine.
    332 !>
     348!> Performs one pass through data as part of multiple real fft (fourier analysis) routine.
     349!>
    333350!> Method:
    334 !> 
     351!>
    335352!> a       is first real input vector
    336353!>         equivalence b(1) with a(ifac*la*inc1+1)
     
    351368!>         2 - ifac not catered for
    352369!>         3 - ifac only catered for if la=n/ifac
    353 !------------------------------------------------------------------------------!
     370!--------------------------------------------------------------------------------------------------!
    354371 SUBROUTINE qpassm( a, b, c, d, trigs, inc1, inc2, inc3, inc4, lot, n, ifac, la, ierr )
    355372
    356373    USE kinds
    357374
    358     IMPLICIT NONE
    359 
    360     INTEGER(iwp) ::  ierr !<
    361     INTEGER(iwp) ::  ifac !<
    362     INTEGER(iwp) ::  inc1 !<
    363     INTEGER(iwp) ::  inc2 !<
    364     INTEGER(iwp) ::  inc3 !<
    365     INTEGER(iwp) ::  inc4 !<
    366     INTEGER(iwp) ::  la   !<
    367     INTEGER(iwp) ::  lot  !<
    368     INTEGER(iwp) ::  n    !<
     375    IMPLICIT NONE
     376
     377    INTEGER(iwp) ::  ierr  !<
     378    INTEGER(iwp) ::  ifac  !<
     379    INTEGER(iwp) ::  inc1  !<
     380    INTEGER(iwp) ::  inc2  !<
     381    INTEGER(iwp) ::  inc3  !<
     382    INTEGER(iwp) ::  inc4  !<
     383    INTEGER(iwp) ::  la    !<
     384    INTEGER(iwp) ::  lot   !<
     385    INTEGER(iwp) ::  n     !<
     386
     387    INTEGER(iwp) ::  i      !<
     388    INTEGER(iwp) ::  ia     !<
     389    INTEGER(iwp) ::  ib     !<
     390    INTEGER(iwp) ::  ibase  !<
     391    INTEGER(iwp) ::  ic     !<
     392    INTEGER(iwp) ::  id     !<
     393    INTEGER(iwp) ::  ie     !<
     394    INTEGER(iwp) ::  if     !<
     395    INTEGER(iwp) ::  ig     !<
     396    INTEGER(iwp) ::  igo    !<
     397    INTEGER(iwp) ::  ih     !<
     398    INTEGER(iwp) ::  iink   !<
     399    INTEGER(iwp) ::  ijk    !<
     400    INTEGER(iwp) ::  ijump  !<
     401    INTEGER(iwp) ::  j      !<
     402    INTEGER(iwp) ::  ja     !<
     403    INTEGER(iwp) ::  jb     !<
     404    INTEGER(iwp) ::  jbase  !<
     405    INTEGER(iwp) ::  jc     !<
     406    INTEGER(iwp) ::  jd     !<
     407    INTEGER(iwp) ::  je     !<
     408    INTEGER(iwp) ::  jf     !<
     409    INTEGER(iwp) ::  jink   !<
     410    INTEGER(iwp) ::  k      !<
     411    INTEGER(iwp) ::  kb     !<
     412    INTEGER(iwp) ::  kc     !<
     413    INTEGER(iwp) ::  kd     !<
     414    INTEGER(iwp) ::  ke     !<
     415    INTEGER(iwp) ::  kf     !<
     416    INTEGER(iwp) ::  kstop  !<
     417    INTEGER(iwp) ::  l      !<
     418    INTEGER(iwp) ::  m      !<
    369419
    370420!
    371421!-- Arrays are dimensioned with n
    372     REAL(wp) ::  a(*)     !<
    373     REAL(wp) ::  b(*)     !<
    374     REAL(wp) ::  c(*)     !<
    375     REAL(wp) ::  d(*)     !<
    376     REAL(wp) ::  trigs(*) !<
    377  
    378     REAL(wp) ::  a0     !<
    379     REAL(wp) ::  a1     !<
    380     REAL(wp) ::  a10    !<
    381     REAL(wp) ::  a11    !<
    382     REAL(wp) ::  a2     !<
    383     REAL(wp) ::  a20    !<
    384     REAL(wp) ::  a21    !<
    385     REAL(wp) ::  a3     !<
    386     REAL(wp) ::  a4     !<
    387     REAL(wp) ::  a5     !<
    388     REAL(wp) ::  a6     !<
    389     REAL(wp) ::  b0     !<
    390     REAL(wp) ::  b1     !<
    391     REAL(wp) ::  b10    !<
    392     REAL(wp) ::  b11    !<
    393     REAL(wp) ::  b2     !<
    394     REAL(wp) ::  b20    !<
    395     REAL(wp) ::  b21    !<
    396     REAL(wp) ::  b3     !<
    397     REAL(wp) ::  b4     !<
    398     REAL(wp) ::  b5     !<
    399     REAL(wp) ::  b6     !<
    400     REAL(wp) ::  c1     !<
    401     REAL(wp) ::  c2     !<
    402     REAL(wp) ::  c3     !<
    403     REAL(wp) ::  c4     !<
    404     REAL(wp) ::  c5     !<
    405     REAL(wp) ::  qrt5   !<
    406     REAL(wp) ::  s1     !<
    407     REAL(wp) ::  s2     !<
    408     REAL(wp) ::  s3     !<
    409     REAL(wp) ::  s4     !<
    410     REAL(wp) ::  s5     !<
    411     REAL(wp) ::  sin36  !<
    412     REAL(wp) ::  sin45  !<
    413     REAL(wp) ::  sin60  !<
    414     REAL(wp) ::  sin72  !<
    415     REAL(wp) ::  z      !<
    416     REAL(wp) ::  zqrt5  !<
    417     REAL(wp) ::  zsin36 !<
    418     REAL(wp) ::  zsin45 !<
    419     REAL(wp) ::  zsin60 !<
    420     REAL(wp) ::  zsin72 !<
    421 
    422     INTEGER(iwp) ::  i     !<
    423     INTEGER(iwp) ::  ia    !<
    424     INTEGER(iwp) ::  ib    !<
    425     INTEGER(iwp) ::  ibase !<
    426     INTEGER(iwp) ::  ic    !<
    427     INTEGER(iwp) ::  id    !<
    428     INTEGER(iwp) ::  ie    !<
    429     INTEGER(iwp) ::  if    !<
    430     INTEGER(iwp) ::  ig    !<
    431     INTEGER(iwp) ::  igo   !<
    432     INTEGER(iwp) ::  ih    !<
    433     INTEGER(iwp) ::  iink  !<
    434     INTEGER(iwp) ::  ijk   !<
    435     INTEGER(iwp) ::  ijump !<
    436     INTEGER(iwp) ::  j     !<
    437     INTEGER(iwp) ::  ja    !<
    438     INTEGER(iwp) ::  jb    !<
    439     INTEGER(iwp) ::  jbase !<
    440     INTEGER(iwp) ::  jc    !<
    441     INTEGER(iwp) ::  jd    !<
    442     INTEGER(iwp) ::  je    !<
    443     INTEGER(iwp) ::  jf    !<
    444     INTEGER(iwp) ::  jink  !<
    445     INTEGER(iwp) ::  k     !<
    446     INTEGER(iwp) ::  kb    !<
    447     INTEGER(iwp) ::  kc    !<
    448     INTEGER(iwp) ::  kd    !<
    449     INTEGER(iwp) ::  ke    !<
    450     INTEGER(iwp) ::  kf    !<
    451     INTEGER(iwp) ::  kstop !<
    452     INTEGER(iwp) ::  l     !<
    453     INTEGER(iwp) ::  m     !<
     422    REAL(wp) ::  a(*)      !<
     423    REAL(wp) ::  b(*)      !<
     424    REAL(wp) ::  c(*)      !<
     425    REAL(wp) ::  d(*)      !<
     426    REAL(wp) ::  trigs(*)  !<
     427
     428    REAL(wp) ::  a0      !<
     429    REAL(wp) ::  a1      !<
     430    REAL(wp) ::  a10     !<
     431    REAL(wp) ::  a11     !<
     432    REAL(wp) ::  a2      !<
     433    REAL(wp) ::  a20     !<
     434    REAL(wp) ::  a21     !<
     435    REAL(wp) ::  a3      !<
     436    REAL(wp) ::  a4      !<
     437    REAL(wp) ::  a5      !<
     438    REAL(wp) ::  a6      !<
     439    REAL(wp) ::  b0      !<
     440    REAL(wp) ::  b1      !<
     441    REAL(wp) ::  b10     !<
     442    REAL(wp) ::  b11     !<
     443    REAL(wp) ::  b2      !<
     444    REAL(wp) ::  b20     !<
     445    REAL(wp) ::  b21     !<
     446    REAL(wp) ::  b3      !<
     447    REAL(wp) ::  b4      !<
     448    REAL(wp) ::  b5      !<
     449    REAL(wp) ::  b6      !<
     450    REAL(wp) ::  c1      !<
     451    REAL(wp) ::  c2      !<
     452    REAL(wp) ::  c3      !<
     453    REAL(wp) ::  c4      !<
     454    REAL(wp) ::  c5      !<
     455    REAL(wp) ::  qrt5    !<
     456    REAL(wp) ::  s1      !<
     457    REAL(wp) ::  s2      !<
     458    REAL(wp) ::  s3      !<
     459    REAL(wp) ::  s4      !<
     460    REAL(wp) ::  s5      !<
     461    REAL(wp) ::  sin36   !<
     462    REAL(wp) ::  sin45   !<
     463    REAL(wp) ::  sin60   !<
     464    REAL(wp) ::  sin72   !<
     465    REAL(wp) ::  z       !<
     466    REAL(wp) ::  zqrt5   !<
     467    REAL(wp) ::  zsin36  !<
     468    REAL(wp) ::  zsin45  !<
     469    REAL(wp) ::  zsin60  !<
     470    REAL(wp) ::  zsin72  !<
    454471
    455472
     
    468485    iink  = la * inc1
    469486    jink  = la * inc2
    470     ijump = (ifac-1) * iink
    471     kstop = ( n-ifac ) / ( 2*ifac )
     487    ijump = (ifac - 1) * iink
     488    kstop = ( n - ifac ) / ( 2 * ifac )
    472489
    473490    ibase = 0
     
    490507          ib = ia + iink
    491508          ja = 1
    492           jb = ja + (2*m-la) * inc2
     509          jb = ja + ( 2 * m - la ) * inc2
    493510
    494511          IF ( la /= m )  THEN
     
    510527             ENDDO
    511528
    512              ja = ja + jink
    513              jink = 2 * jink
    514              jb = jb - jink
     529             ja    = ja + jink
     530             jink  = 2 * jink
     531             jb    = jb - jink
    515532             ibase = ibase + ijump
    516533             ijump = 2 * ijump + iink
     
    530547                      !OCL NOVREC
    531548                      DO  ijk = 1, lot
    532                          c(ja+j) = a(ia+i) + (c1*a(ib+i)+s1*b(ib+i))
    533                          c(jb+j) = a(ia+i) - (c1*a(ib+i)+s1*b(ib+i))
    534                          d(ja+j) = (c1*b(ib+i)-s1*a(ib+i)) + b(ia+i)
    535                          d(jb+j) = (c1*b(ib+i)-s1*a(ib+i)) - b(ia+i)
     549                         c(ja+j) = a(ia+i) + ( c1 * a(ib+i) + s1 * b(ib+i) )
     550                         c(jb+j) = a(ia+i) - ( c1 * a(ib+i) + s1 * b(ib+i) )
     551                         d(ja+j) = ( c1 * b(ib+i) - s1 * a(ib+i) ) + b(ia+i)
     552                         d(jb+j) = ( c1 * b(ib+i) - s1 * a(ib+i) ) - b(ia+i)
    536553                         i = i + inc3
    537554                         j = j + inc4
     
    562579                   i = i + inc3
    563580                   j = j + inc4
    564 
    565581                ENDDO
    566582                ibase = ibase + inc1
     
    570586          ELSE
    571587
    572              z = 1.0_wp/REAL(n,KIND=wp)
     588             z = 1.0_wp / REAL( n, KIND = wp )
    573589             DO  l = 1, la
    574590                i = ibase
     
    578594                !OCL NOVREC
    579595                DO  ijk = 1, lot
    580                    c(ja+j) = z*(a(ia+i)+a(ib+i))
    581                    c(jb+j) = z*(a(ia+i)-a(ib+i))
     596                   c(ja+j) = z * ( a(ia+i) + a(ib+i) )
     597                   c(jb+j) = z * ( a(ia+i) - a(ib+i) )
    582598                   i = i + inc3
    583599                   j = j + inc4
     
    597613          ic = ib + iink
    598614          ja = 1
    599           jb = ja + (2*m-la) * inc2
     615          jb = ja + ( 2 * m - la) * inc2
    600616          jc = jb
    601617
     
    609625                !OCL NOVREC
    610626                DO  ijk = 1, lot
    611                    c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
    612                    c(jb+j) = a(ia+i) - 0.5_wp*(a(ib+i)+a(ic+i))
    613                    d(jb+j) = sin60*(a(ic+i)-a(ib+i))
     627                   c(ja+j) = a(ia+i) + ( a(ib+i) + a(ic+i) )
     628                   c(jb+j) = a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) )
     629                   d(jb+j) = sin60 * ( a(ic+i) - a(ib+i) )
    614630                   i = i + inc3
    615631                   j = j + inc4
     
    644660                      !OCL NOVREC
    645661                      DO  ijk = 1, lot
    646                          a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c2*a(ic+i)+s2*b(ic+i))
    647                          b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c2*b(ic+i)-s2*a(ic+i))
    648                          a2 = a(ia+i) - 0.5_wp*a1
    649                          b2 = b(ia+i) - 0.5_wp*b1
    650                          a3 = sin60*((c1*a(ib+i)+s1*b(ib+i))-(c2*a(ic+i)+s2*b(ic+i)))
    651                          b3 = sin60*((c1*b(ib+i)-s1*a(ib+i))-(c2*b(ic+i)-s2*a(ic+i)))
     662                         a1 = ( c1 * a(ib+i) + s1 * b(ib+i) ) + ( c2 * a(ic+i) + s2 * b(ic+i) )
     663                         b1 = ( c1 * b(ib+i) - s1 * a(ib+i) ) + ( c2 * b(ic+i) - s2 * a(ic+i) )
     664                         a2 = a(ia+i) - 0.5_wp * a1
     665                         b2 = b(ia+i) - 0.5_wp * b1
     666                         a3 = sin60 * ( ( c1 * a(ib+i) + s1 * b(ib+i) ) - ( c2 * a(ic+i) + s2      &
     667                              * b(ic+i) ) )
     668                         b3 = sin60 * ( ( c1 * b(ib+i) - s1 * a(ib+i) ) - ( c2 * b(ic+i) - s2      &
     669                              * a(ic+i) ) )
    652670                         c(ja+j) = a(ia+i) + a1
    653671                         d(ja+j) = b(ia+i) + b1
     
    655673                         d(jb+j) = b2 - a3
    656674                         c(jc+j) = a2 - b3
    657                          d(jc+j) = -(b2+a3)
     675                         d(jc+j) = - ( b2 + a3 )
    658676                         i = i + inc3
    659677                         j = j + inc4
     
    680698                !OCL NOVREC
    681699                DO  ijk = 1, lot
    682                    c(ja+j) = a(ia+i) + 0.5_wp*(a(ib+i)-a(ic+i))
    683                    d(ja+j) = -sin60*(a(ib+i)+a(ic+i))
    684                    c(jb+j) = a(ia+i) - (a(ib+i)-a(ic+i))
     700                   c(ja+j) = a(ia+i) + 0.5_wp * ( a(ib+i) - a(ic+i) )
     701                   d(ja+j) = -sin60 * ( a(ib+i) + a(ic+i) )
     702                   c(jb+j) = a(ia+i) - ( a(ib+i) - a(ic+i) )
    685703                   i = i + inc3
    686704                   j = j + inc4
     
    692710          ELSE
    693711
    694              z = 1.0_wp / REAL( n, KIND=wp )
    695              zsin60 = z*sin60
     712             z = 1.0_wp / REAL( n, KIND = wp )
     713             zsin60 = z * sin60
    696714             DO  l = 1, la
    697715                i = ibase
     
    701719                !OCL NOVREC
    702720                DO  ijk = 1, lot
    703                    c(ja+j) = z*(a(ia+i)+(a(ib+i)+a(ic+i)))
    704                    c(jb+j) = z*(a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))
    705                    d(jb+j) = zsin60*(a(ic+i)-a(ib+i))
     721                   c(ja+j) = z * ( a(ia+i) + ( a (ib+i) + a(ic+i) ) )
     722                   c(jb+j) = z * ( a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) ) )
     723                   d(jb+j) = zsin60 * ( a(ic+i) -a(ib+i) )
    706724                   i = i + inc3
    707725                   j = j + inc4
     
    722740          id = ic + iink
    723741          ja = 1
    724           jb = ja + (2*m-la) * inc2
    725           jc = jb + 2*m*inc2
     742          jb = ja + ( 2 * m - la ) * inc2
     743          jc = jb + 2 * m * inc2
    726744          jd = jb
    727745
     
    736754                !OCL NOVREC
    737755                DO  ijk = 1, lot
    738                    c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))
    739                    c(jc+j) = (a(ia+i)+a(ic+i)) - (a(ib+i)+a(id+i))
     756                   c(ja+j) = ( a(ia+i) + a(ic+i) ) + ( a(ib+i) + a(id+i) )
     757                   c(jc+j) = ( a(ia+i) + a(ic+i) ) - ( a(ib+i) + a(id+i) )
    740758                   c(jb+j) = a(ia+i) - a(ic+i)
    741759                   d(jb+j) = a(id+i) - a(ib+i)
     
    747765             ENDDO
    748766
    749              ja = ja + jink
    750              jink = 2 * jink
    751              jb = jb + jink
    752              jc = jc - jink
    753              jd = jd - jink
     767             ja    = ja + jink
     768             jink  = 2 * jink
     769             jb    = jb + jink
     770             jc    = jc - jink
     771             jd    = jd - jink
    754772             ibase = ibase + ijump
    755773             ijump = 2 * ijump + iink
     
    775793                      !OCL NOVREC
    776794                      DO  ijk = 1, lot
    777                          a0 = a(ia+i) + (c2*a(ic+i)+s2*b(ic+i))
    778                          a2 = a(ia+i) - (c2*a(ic+i)+s2*b(ic+i))
    779                          a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c3*a(id+i)+s3*b(id+i))
    780                          a3 = (c1*a(ib+i)+s1*b(ib+i)) - (c3*a(id+i)+s3*b(id+i))
    781                          b0 = b(ia+i) + (c2*b(ic+i)-s2*a(ic+i))
    782                          b2 = b(ia+i) - (c2*b(ic+i)-s2*a(ic+i))
    783                          b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c3*b(id+i)-s3*a(id+i))
    784                          b3 = (c1*b(ib+i)-s1*a(ib+i)) - (c3*b(id+i)-s3*a(id+i))
     795                         a0 = a(ia+i) + ( c2 * a(ic+i) + s2 * b(ic+i) )
     796                         a2 = a(ia+i) - ( c2 * a(ic+i) + s2 * b(ic+i) )
     797                         a1 = ( c1 * a(ib+i) + s1 * b(ib+i) ) + ( c3 * a(id+i) + s3 * b(id+i) )
     798                         a3 = ( c1 * a(ib+i) + s1 * b(ib+i) ) - ( c3 * a(id+i) + s3 * b(id+i) )
     799                         b0 = b(ia+i) + ( c2 * b(ic+i) - s2 * a(ic+i) )
     800                         b2 = b(ia+i) - ( c2 * b(ic+i) - s2 * a(ic+i) )
     801                         b1 = ( c1 * b(ib+i) - s1 * a(ib+i) ) + ( c3 * b(id+i) - s3 * a(id+i) )
     802                         b3 = ( c1 * b(ib+i) - s1 * a(ib+i) ) - ( c3 * b(id+i) - s3 * a(id+i) )
    785803                         c(ja+j) = a0 + a1
    786804                         c(jc+j) = a0 - a1
     
    790808                         c(jd+j) = a2 - b3
    791809                         d(jb+j) = b2 - a3
    792                          d(jd+j) = -(b2+a3)
     810                         d(jd+j) = - ( b2 + a3 )
    793811                         i = i + inc3
    794812                         j = j + inc4
     
    817835                !OCL NOVREC
    818836                DO  ijk = 1, lot
    819                    c(ja+j) = a(ia+i) + sin45*(a(ib+i)-a(id+i))
    820                    c(jb+j) = a(ia+i) - sin45*(a(ib+i)-a(id+i))
    821                    d(ja+j) = -a(ic+i) - sin45*(a(ib+i)+a(id+i))
    822                    d(jb+j) = a(ic+i) - sin45*(a(ib+i)+a(id+i))
     837                   c(ja+j) = a(ia+i) + sin45 * ( a(ib+i) - a(id+i) )
     838                   c(jb+j) = a(ia+i) - sin45 * ( a(ib+i) - a(id+i) )
     839                   d(ja+j) = - a(ic+i) - sin45 * ( a(ib+i) + a(id+i) )
     840                   d(jb+j) = a(ic+i) - sin45 * ( a(ib+i) + a(id+i) )
    823841                   i = i + inc3
    824842                   j = j + inc4
     
    830848          ELSE
    831849
    832              z = 1.0_wp / REAL( n, KIND=wp )
     850             z = 1.0_wp / REAL( n, KIND = wp )
    833851             DO  l = 1, la
    834852                i = ibase
     
    838856                !OCL NOVREC
    839857                DO  ijk = 1, lot
    840                    c(ja+j) = z*((a(ia+i)+a(ic+i))+(a(ib+i)+a(id+i)))
    841                    c(jc+j) = z*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i)))
    842                    c(jb+j) = z*(a(ia+i)-a(ic+i))
    843                    d(jb+j) = z*(a(id+i)-a(ib+i))
     858                   c(ja+j) = z * ( ( a(ia+i) + a(ic+i) ) + ( a(ib+i) + a(id+i) ) )
     859                   c(jc+j) = z * ( ( a(ia+i) + a(ic+i) ) - (a (ib+i) + a(id+i) ) )
     860                   c(jb+j) = z * ( a(ia+i) - a(ic+i) )
     861                   d(jb+j) = z * ( a(id+i) - a(ib+i) )
    844862                   i = i + inc3
    845863                   j = j + inc4
     
    861879          ie = id + iink
    862880          ja = 1
    863           jb = ja + (2*m-la) * inc2
    864           jc = jb + 2*m*inc2
     881          jb = ja + ( 2 * m - la ) * inc2
     882          jc = jb + 2 * m * inc2
    865883          jd = jc
    866884          je = jb
     
    879897                   a2 = a(ic+i) + a(id+i)
    880898                   a4 = a(ic+i) - a(id+i)
    881                    a5 = a(ia+i) - 0.25_wp*(a1+a2)
    882                    a6 = qrt5*(a1-a2)
    883                    c(ja+j) = a(ia+i) + (a1+a2)
     899                   a5 = a(ia+i) - 0.25_wp * ( a1 + a2 )
     900                   a6 = qrt5 * ( a1 - a2 )
     901                   c(ja+j) = a(ia+i) + ( a1 + a2 )
    884902                   c(jb+j) = a5 + a6
    885903                   c(jc+j) = a5 - a6
    886                    d(jb+j) = -sin72*a3 - sin36*a4
    887                    d(jc+j) = -sin36*a3 + sin72*a4
     904                   d(jb+j) = - sin72 * a3 - sin36 * a4
     905                   d(jc+j) = - sin36 * a3 + sin72 * a4
    888906                   i = i + inc3
    889907                   j = j + inc4
     
    925943                      !OCL NOVREC
    926944                      DO  ijk = 1, lot
    927                          a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c4*a(ie+i)+s4*b(ie+i))
    928                          a3 = (c1*a(ib+i)+s1*b(ib+i)) - (c4*a(ie+i)+s4*b(ie+i))
    929                          a2 = (c2*a(ic+i)+s2*b(ic+i)) + (c3*a(id+i)+s3*b(id+i))
    930                          a4 = (c2*a(ic+i)+s2*b(ic+i)) - (c3*a(id+i)+s3*b(id+i))
    931                          b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c4*b(ie+i)-s4*a(ie+i))
    932                          b3 = (c1*b(ib+i)-s1*a(ib+i)) - (c4*b(ie+i)-s4*a(ie+i))
    933                          b2 = (c2*b(ic+i)-s2*a(ic+i)) + (c3*b(id+i)-s3*a(id+i))
    934                          b4 = (c2*b(ic+i)-s2*a(ic+i)) - (c3*b(id+i)-s3*a(id+i))
    935                          a5 = a(ia+i) - 0.25_wp*(a1+a2)
    936                          a6 = qrt5*(a1-a2)
    937                          b5 = b(ia+i) - 0.25_wp*(b1+b2)
    938                          b6 = qrt5*(b1-b2)
     945                         a1  = (c1 * a(ib+i) + s1 * b(ib+i) ) + ( c4 * a(ie+i) + s4 * b(ie+i) )
     946                         a3  = (c1 * a(ib+i) + s1 * b(ib+i) ) - ( c4 * a(ie+i) + s4 * b(ie+i) )
     947                         a2  = (c2 * a(ic+i) + s2 * b(ic+i) ) + ( c3 * a(id+i) + s3 * b(id+i) )
     948                         a4  = (c2 * a(ic+i) + s2 * b(ic+i) ) - ( c3 * a(id+i) + s3 * b(id+i) )
     949                         b1  = (c1 * b(ib+i) - s1 * a(ib+i) ) + ( c4 * b(ie+i) - s4 * a(ie+i) )
     950                         b3  = (c1 * b(ib+i) - s1 * a(ib+i) ) - ( c4 * b(ie+i) - s4 * a(ie+i) )
     951                         b2  = (c2 * b(ic+i) - s2 * a(ic+i) ) + ( c3 * b(id+i) - s3 * a(id+i) )
     952                         b4  = (c2 * b(ic+i) - s2 * a(ic+i) ) - ( c3 * b(id+i) - s3 * a(id+i) )
     953                         a5  = a(ia+i) - 0.25_wp * ( a1 + a2 )
     954                         a6  = qrt5 * ( a1 - a2 )
     955                         b5  = b(ia+i) - 0.25_wp * ( b1 + b2 )
     956                         b6  = qrt5 * ( b1 - b2 )
    939957                         a10 = a5 + a6
    940958                         a20 = a5 - a6
    941959                         b10 = b5 + b6
    942960                         b20 = b5 - b6
    943                          a11 = sin72*b3 + sin36*b4
    944                          a21 = sin36*b3 - sin72*b4
    945                          b11 = sin72*a3 + sin36*a4
    946                          b21 = sin36*a3 - sin72*a4
    947                          c(ja+j) = a(ia+i) + (a1+a2)
     961                         a11 = sin72 * b3 + sin36 * b4
     962                         a21 = sin36 * b3 - sin72 * b4
     963                         b11 = sin72 * a3 + sin36 * a4
     964                         b21 = sin36 * a3 - sin72 * a4
     965                         c(ja+j) = a(ia+i) + ( a1 + a2 )
    948966                         c(jb+j) = a10 + a11
    949967                         c(je+j) = a10 - a11
    950968                         c(jc+j) = a20 + a21
    951969                         c(jd+j) = a20 - a21
    952                          d(ja+j) = b(ia+i) + (b1+b2)
     970                         d(ja+j) = b(ia+i) + ( b1 + b2 )
    953971                         d(jb+j) = b10 - b11
    954                          d(je+j) = -(b10+b11)
     972                         d(je+j) = - ( b10 + b11 )
    955973                         d(jc+j) = b20 - b21
    956                          d(jd+j) = -(b20+b21)
     974                         d(jd+j) = - ( b20 + b21 )
    957975                         i = i + inc3
    958976                         j = j + inc4
     
    9851003                   a2 = a(ic+i) + a(id+i)
    9861004                   a4 = a(ic+i) - a(id+i)
    987                    a5 = a(ia+i) + 0.25_wp*(a3-a4)
    988                    a6 = qrt5*(a3+a4)
     1005                   a5 = a(ia+i) + 0.25_wp * ( a3 - a4 )
     1006                   a6 = qrt5 * ( a3 + a4 )
    9891007                   c(ja+j) = a5 + a6
    9901008                   c(jb+j) = a5 - a6
    991                    c(jc+j) = a(ia+i) - (a3-a4)
    992                    d(ja+j) = -sin36*a1 - sin72*a2
    993                    d(jb+j) = -sin72*a1 + sin36*a2
     1009                   c(jc+j) = a(ia+i) - ( a3 - a4 )
     1010                   d(ja+j) = - sin36 * a1 - sin72 * a2
     1011                   d(jb+j) = - sin72 * a1 + sin36 * a2
    9941012                   i = i + inc3
    9951013                   j = j + inc4
     
    10011019          ELSE
    10021020
    1003              z = 1.0_wp / REAL( n, KIND=wp )
     1021             z = 1.0_wp / REAL( n, KIND = wp )
    10041022             zqrt5  = z * qrt5
    10051023             zsin36 = z * sin36
     
    10161034                   a2 = a(ic+i) + a(id+i)
    10171035                   a4 = a(ic+i) - a(id+i)
    1018                    a5 = z*(a(ia+i)-0.25_wp*(a1+a2))
    1019                    a6 = zqrt5*(a1-a2)
    1020                    c(ja+j) = z*(a(ia+i)+(a1+a2))
     1036                   a5 = z * ( a (ia+i) - 0.25_wp * ( a1 + a2 ) )
     1037                   a6 = zqrt5 * ( a1 - a2 )
     1038                   c(ja+j) = z * ( a(ia+i) + (a1+a2) )
    10211039                   c(jb+j) = a5 + a6
    10221040                   c(jc+j) = a5 - a6
    1023                    d(jb+j) = -zsin72*a3 - zsin36*a4
    1024                    d(jc+j) = -zsin36*a3 + zsin72*a4
     1041                   d(jb+j) = - zsin72 * a3 - zsin36 * a4
     1042                   d(jc+j) = - zsin36 * a3 + zsin72 * a4
    10251043                   i = i + inc3
    10261044                   j = j + inc4
     
    10431061          if = ie + iink
    10441062          ja = 1
    1045           jb = ja + (2*m-la) * inc2
    1046           jc = jb + 2*m*inc2
    1047           jd = jc + 2*m*inc2
     1063          jb = ja + ( 2 * m - la ) * inc2
     1064          jc = jb + 2 * m * inc2
     1065          jd = jc + 2 * m * inc2
    10481066          je = jc
    10491067          jf = jb
     
    10581076                !OCL NOVREC
    10591077                DO  ijk = 1, lot
    1060                    a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
    1061                    c(ja+j) = (a(ia+i)+a(id+i)) + a11
    1062                    c(jc+j) = (a(ia+i)+a(id+i)-0.5_wp*a11)
    1063                    d(jc+j) = sin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
    1064                    a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
    1065                    c(jb+j) = (a(ia+i)-a(id+i)) - 0.5_wp*a11
    1066                    d(jb+j) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
    1067                    c(jd+j) = (a(ia+i)-a(id+i)) + a11
     1078                   a11 = ( a (ic+i) + a(if+i) ) + ( a(ib+i) + a(ie+i) )
     1079                   c(ja+j) = ( a(ia+i) + a(id+i) ) + a11
     1080                   c(jc+j) = ( a(ia+i) + a(id+i) - 0.5_wp * a11 )
     1081                   d(jc+j) = sin60 * ( ( a(ic+i) + a(if+i) ) - ( a(ib+i) + a(ie+i) ) )
     1082                   a11 = ( a (ic+i) - a(if+i) ) + ( a(ie+i) - a(ib+i) )
     1083                   c(jb+j) = ( a(ia+i) - a(id+i) ) - 0.5_wp * a11
     1084                   d(jb+j) = sin60 * ( ( a(ie+i) - a(ib+i) ) - ( a(ic+i) - a(if+i) ) )
     1085                   c(jd+j) = ( a (ia+i) - a(id+i) ) + a11
    10681086                   i = i + inc3
    10691087                   j = j + inc4
     
    10721090                jbase = jbase + inc2
    10731091             ENDDO
    1074              ja = ja + jink
    1075              jink = 2 * jink
    1076              jb = jb + jink
    1077              jc = jc + jink
    1078              jd = jd - jink
    1079              je = je - jink
    1080              jf = jf - jink
     1092             ja    = ja + jink
     1093             jink  = 2 * jink
     1094             jb    = jb + jink
     1095             jc    = jc + jink
     1096             jd    = jd - jink
     1097             je    = je - jink
     1098             jf    = jf - jink
    10811099             ibase = ibase + ijump
    10821100             ijump = 2 * ijump + iink
     
    11081126                      !OCL NOVREC
    11091127                      DO  ijk = 1, lot
    1110                          a1 = c1*a(ib+i) + s1*b(ib+i)
    1111                          b1 = c1*b(ib+i) - s1*a(ib+i)
    1112                          a2 = c2*a(ic+i) + s2*b(ic+i)
    1113                          b2 = c2*b(ic+i) - s2*a(ic+i)
    1114                          a3 = c3*a(id+i) + s3*b(id+i)
    1115                          b3 = c3*b(id+i) - s3*a(id+i)
    1116                          a4 = c4*a(ie+i) + s4*b(ie+i)
    1117                          b4 = c4*b(ie+i) - s4*a(ie+i)
    1118                          a5 = c5*a(if+i) + s5*b(if+i)
    1119                          b5 = c5*b(if+i) - s5*a(if+i)
    1120                          a11 = (a2+a5) + (a1+a4)
    1121                          a20 = (a(ia+i)+a3) - 0.5_wp*a11
    1122                          a21 = sin60*((a2+a5)-(a1+a4))
    1123                          b11 = (b2+b5) + (b1+b4)
    1124                          b20 = (b(ia+i)+b3) - 0.5_wp*b11
    1125                          b21 = sin60*((b2+b5)-(b1+b4))
    1126                          c(ja+j) = (a(ia+i)+a3) + a11
    1127                          d(ja+j) = (b(ia+i)+b3) + b11
     1128                         a1  = c1 * a(ib+i) + s1 * b(ib+i)
     1129                         b1  = c1 * b(ib+i) - s1 * a(ib+i)
     1130                         a2  = c2 * a(ic+i) + s2 * b(ic+i)
     1131                         b2  = c2 * b(ic+i) - s2 * a(ic+i)
     1132                         a3  = c3 * a(id+i) + s3 * b(id+i)
     1133                         b3  = c3 * b(id+i) - s3 * a(id+i)
     1134                         a4  = c4 * a(ie+i) + s4 * b(ie+i)
     1135                         b4  = c4 * b(ie+i) - s4 * a(ie+i)
     1136                         a5  = c5 * a(if+i) + s5 * b(if+i)
     1137                         b5  = c5 * b(if+i) - s5 * a(if+i)
     1138                         a11 = ( a2 + a5 ) + ( a1 + a4 )
     1139                         a20 = ( a(ia+i) + a3 ) - 0.5_wp * a11
     1140                         a21 = sin60 * ( ( a2 + a5 ) - ( a1 + a4 ) )
     1141                         b11 = ( b2 + b5 ) + ( b1 + b4 )
     1142                         b20 = ( b(ia+i) + b3 ) - 0.5_wp * b11
     1143                         b21 = sin60 * ( ( b2 + b5 ) - ( b1 + b4 ) )
     1144                         c(ja+j) = ( a(ia+i) + a3 ) + a11
     1145                         d(ja+j) = ( b(ia+i) + b3 ) + b11
    11281146                         c(jc+j) = a20 - b21
    11291147                         d(jc+j) = a21 + b20
    11301148                         c(je+j) = a20 + b21
    11311149                         d(je+j) = a21 - b20
    1132                          a11 = (a2-a5) + (a4-a1)
    1133                          a20 = (a(ia+i)-a3) - 0.5_wp*a11
    1134                          a21 = sin60*((a4-a1)-(a2-a5))
    1135                          b11 = (b5-b2) - (b4-b1)
    1136                          b20 = (b3-b(ia+i)) - 0.5_wp*b11
    1137                          b21 = sin60*((b5-b2)+(b4-b1))
     1150                         a11 = ( a2 - a5 ) + ( a4 - a1 )
     1151                         a20 = ( a(ia+i) - a3 ) - 0.5_wp * a11
     1152                         a21 = sin60 * ( ( a4 - a1 ) - ( a2 - a5 ) )
     1153                         b11 = ( b5 - b2 ) - ( b4 - b1 )
     1154                         b20 = ( b3 - b(ia+i) ) - 0.5_wp * b11
     1155                         b21 = sin60 * ( ( b5 - b2 ) + ( b4 - b1 ) )
    11381156                         c(jb+j) = a20 - b21
    11391157                         d(jb+j) = a21 - b20
    1140                          c(jd+j) = a11 + (a(ia+i)-a3)
    1141                          d(jd+j) = b11 + (b3-b(ia+i))
     1158                         c(jd+j) = a11 + ( a(ia+i) - a3 )
     1159                         d(jd+j) = b11 + ( b3 - b(ia+i) )
    11421160                         c(jf+j) = a20 + b21
    11431161                         d(jf+j) = a21 + b20
     
    11691187                !OCL NOVREC
    11701188                DO  ijk = 1, lot
    1171                    c(ja+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) + sin60*(a(ib+i)-a(if+i))
    1172                    d(ja+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) - sin60*(a(ic+i)+a(ie+i))
    1173                    c(jb+j) = a(ia+i) - (a(ic+i)-a(ie+i))
    1174                    d(jb+j) = a(id+i) - (a(ib+i)+a(if+i))
    1175                    c(jc+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) - sin60*(a(ib+i)-a(if+i))
    1176                    d(jc+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) + sin60*(a(ic+i)+a(ie+i))
     1189                   c(ja+j) = ( a (ia+i) + 0.5_wp * ( a(ic+i) - a(ie+i) ) ) + sin60 *               &
     1190                             ( a(ib+i) - a(if+i) )
     1191                   d(ja+j) = - ( a(id+i) + 0.5_wp * ( a(ib+i) + a(if+i) ) ) - sin60 * ( a(ic+i) +  &
     1192                             a(ie+i) )
     1193                   c(jb+j) = a(ia+i) - ( a(ic+i) - a(ie+i) )
     1194                   d(jb+j) = a(id+i) - ( a(ib+i) + a(if+i) )
     1195                   c(jc+j) = ( a(ia+i) + 0.5_wp * ( a(ic+i) - a(ie+i) ) ) - sin60 * ( a(ib+i) -    &
     1196                             a(if+i) )
     1197                   d(jc+j) = - ( a(id+i) + 0.5_wp * ( a(ib+i) + a(if+i) ) ) + sin60 * ( a(ic+i) +  &
     1198                             a(ie+i) )
    11771199                   i = i + inc3
    11781200                   j = j + inc4
     
    11841206          ELSE
    11851207
    1186              z = 1.0_wp/REAL(n,KIND=wp)
    1187              zsin60 = z*sin60
     1208             z = 1.0_wp / REAL( n, KIND = wp )
     1209             zsin60 = z * sin60
    11881210             DO  l = 1, la
    11891211                i = ibase
     
    11931215                !OCL NOVREC
    11941216                DO  ijk = 1, lot
    1195                    a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
    1196                    c(ja+j) = z*((a(ia+i)+a(id+i))+a11)
    1197                    c(jc+j) = z*((a(ia+i)+a(id+i))-0.5_wp*a11)
    1198                    d(jc+j) = zsin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
    1199                    a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
    1200                    c(jb+j) = z*((a(ia+i)-a(id+i))-0.5_wp*a11)
    1201                    d(jb+j) = zsin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
    1202                    c(jd+j) = z*((a(ia+i)-a(id+i))+a11)
     1217                   a11 = ( a(ic+i) + a(if+i) ) + ( a(ib+i) + a(ie+i) )
     1218                   c(ja+j) = z * ( ( a(ia+i) + a(id+i) ) + a11 )
     1219                   c(jc+j) = z * ( ( a(ia+i) + a(id+i) ) - 0.5_wp * a11 )
     1220                   d(jc+j) = zsin60 * ( ( a(ic+i) + a(if+i) ) - ( a(ib+i) + a(ie+i) ) )
     1221                   a11 = ( a(ic+i) - a(if+i) ) + ( a(ie+i) - a(ib+i) )
     1222                   c(jb+j) = z * ( ( a(ia+i) - a(id+i) ) - 0.5_wp * a11 )
     1223                   d(jb+j) = zsin60 * ( ( a(ie+i) - a(ib+i) ) - ( a(ic+i) - a(if+i) ) )
     1224                   c(jd+j) = z * ( ( a(ia+i) - a(id+i) ) + a11 )
    12031225                   i = i + inc3
    12041226                   j = j + inc4
     
    12291251          ja = 1
    12301252          jb = ja + la * inc2
    1231           jc = jb + 2*m*inc2
    1232           jd = jc + 2*m*inc2
    1233           je = jd + 2*m*inc2
    1234           z = 1.0_wp / REAL( n, KIND=wp )
     1253          jc = jb + 2 * m * inc2
     1254          jd = jc + 2 * m * inc2
     1255          je = jd + 2 * m * inc2
     1256          z = 1.0_wp / REAL( n, KIND = wp )
    12351257          zsin45 = z * SQRT( 0.5_wp )
    12361258
     
    12431265             !OCL NOVREC
    12441266             DO  ijk = 1, lot
    1245                 c(ja+j) = z*(((a(ia+i)+a(ie+i))+(a(ic+i)+a(ig+i)))+((a(id+i)+ a(ih+i))+(a(ib+i)+a(if+i))))
    1246                 c(je+j) = z*(((a(ia+i)+a(ie+i))+(a(ic+i)+a(ig+i)))-((a(id+i)+ a(ih+i))+(a(ib+i)+a(if+i))))
    1247                 c(jc+j) = z*((a(ia+i)+a(ie+i))-(a(ic+i)+a(ig+i)))
    1248                 d(jc+j) = z*((a(id+i)+a(ih+i))-(a(ib+i)+a(if+i)))
    1249                 c(jb+j) = z*(a(ia+i)-a(ie+i)) + zsin45*((a(ih+i)-a(id+i))-(a(if+i)-a(ib+i)))
    1250                 c(jd+j) = z*(a(ia+i)-a(ie+i)) - zsin45*((a(ih+i)-a(id+i))-(a(if+i)-a(ib+i)))
    1251                 d(jb+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) + z*(a(ig+i)-a(ic+i))
    1252                 d(jd+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) - z*(a(ig+i)-a(ic+i))
     1267                c(ja+j) = z * ( ( ( a(ia+i) + a(ie+i) ) + ( a(ic+i) + a(ig+i) ) ) + ( ( a(id+i)    &
     1268                          + a(ih+i) ) + ( a(ib+i) + a(if+i) ) ) )
     1269                c(je+j) = z * ( ( ( a(ia+i) + a(ie+i) ) + (a(ic+i) + a(ig+i) ) ) - ( ( a(id+i)     &
     1270                          + a(ih+i) ) + ( a(ib+i) + a(if+i) ) ) )
     1271                c(jc+j) = z * ( ( a(ia+i) + a(ie+i) ) - ( a(ic+i) + a(ig+i) ) )
     1272                d(jc+j) = z * ( ( a(id+i) + a(ih+i) ) - ( a(ib+i) + a(if+i) ) )
     1273                c(jb+j) = z * ( a(ia+i) - a(ie+i) ) + zsin45 * ( ( a(ih+i) - a(id+i) ) - ( a(if+i) &
     1274                          - a(ib+i) ) )
     1275                c(jd+j) = z * ( a(ia+i) - a(ie+i) ) - zsin45 * ( ( a(ih+i) - a(id+i) ) - ( a(if+i) &
     1276                          - a(ib+i) ) )
     1277                d(jb+j) = zsin45 * ( ( a(ih+i) - a(id+i) ) + ( a(if+i) - a(ib+i) ) ) + z *         &
     1278                          ( a(ig+i) - a(ic+i) )
     1279                d(jd+j) = zsin45 * ( ( a(ih+i) - a(id+i) ) + ( a(if+i) - a(ib+i) ) ) - z *         &
     1280                          ( a(ig+i) - a(ic+i) )
    12531281                i = i + inc3
    12541282                j = j + inc4
     
    12621290 END SUBROUTINE qpassm
    12631291
    1264 !------------------------------------------------------------------------------!
     1292!--------------------------------------------------------------------------------------------------!
    12651293! Description:
    12661294! ------------
    12671295!> Same as qpassm, but for backward fft
    1268 !------------------------------------------------------------------------------!
     1296!--------------------------------------------------------------------------------------------------!
    12691297 SUBROUTINE rpassm( a, b, c, d, trigs, inc1, inc2, inc3, inc4, lot, n, ifac, la, ierr )
    12701298
     
    12731301    IMPLICIT NONE
    12741302
    1275     INTEGER(iwp) ::  ierr !<
    1276     INTEGER(iwp) ::  ifac !<
    1277     INTEGER(iwp) ::  inc1 !<
    1278     INTEGER(iwp) ::  inc2 !<
    1279     INTEGER(iwp) ::  inc3 !<
    1280     INTEGER(iwp) ::  inc4 !<
    1281     INTEGER(iwp) ::  la   !<
    1282     INTEGER(iwp) ::  lot  !<
    1283     INTEGER(iwp) ::  n    !<
     1303    INTEGER(iwp) ::  ierr  !<
     1304    INTEGER(iwp) ::  ifac  !<
     1305    INTEGER(iwp) ::  inc1  !<
     1306    INTEGER(iwp) ::  inc2  !<
     1307    INTEGER(iwp) ::  inc3  !<
     1308    INTEGER(iwp) ::  inc4  !<
     1309    INTEGER(iwp) ::  la    !<
     1310    INTEGER(iwp) ::  lot   !<
     1311    INTEGER(iwp) ::  n     !<
    12841312!
    12851313!-- Arrays are dimensioned with n
    1286     REAL(wp) ::  a(*)     !<
    1287     REAL(wp) ::  b(*)     !<
    1288     REAL(wp) ::  c(*)     !<
    1289     REAL(wp) ::  d(*)     !<
    1290     REAL(wp) ::  trigs(*) !<
    1291 
    1292     REAL(wp) ::  c1     !<
    1293     REAL(wp) ::  c2     !<
    1294     REAL(wp) ::  c3     !<
    1295     REAL(wp) ::  c4     !<
    1296     REAL(wp) ::  c5     !<
    1297     REAL(wp) ::  qqrt5  !<
    1298     REAL(wp) ::  qrt5   !<
    1299     REAL(wp) ::  s1     !<
    1300     REAL(wp) ::  s2     !<
    1301     REAL(wp) ::  s3     !<
    1302     REAL(wp) ::  s4     !<
    1303     REAL(wp) ::  s5     !<
    1304     REAL(wp) ::  sin36  !<
    1305     REAL(wp) ::  sin45  !<
    1306     REAL(wp) ::  sin60  !<
    1307     REAL(wp) ::  sin72  !<
    1308     REAL(wp) ::  ssin36 !<
    1309     REAL(wp) ::  ssin45 !<
    1310     REAL(wp) ::  ssin60 !<
    1311     REAL(wp) ::  ssin72 !<
    1312 
    1313     INTEGER(iwp) ::  i     !<
    1314     INTEGER(iwp) ::  ia    !<
    1315     INTEGER(iwp) ::  ib    !<
    1316     INTEGER(iwp) ::  ibase !<
    1317     INTEGER(iwp) ::  ic    !<
    1318     INTEGER(iwp) ::  id    !<
    1319     INTEGER(iwp) ::  ie    !<
    1320     INTEGER(iwp) ::  if    !<
    1321     INTEGER(iwp) ::  igo   !<
    1322     INTEGER(iwp) ::  iink  !<
    1323     INTEGER(iwp) ::  ijk   !<
    1324     INTEGER(iwp) ::  j     !<
    1325     INTEGER(iwp) ::  ja    !<
    1326     INTEGER(iwp) ::  jb    !<
    1327     INTEGER(iwp) ::  jbase !<
    1328     INTEGER(iwp) ::  jc    !<
    1329     INTEGER(iwp) ::  jd    !<
    1330     INTEGER(iwp) ::  je    !<
    1331     INTEGER(iwp) ::  jf    !<
    1332     INTEGER(iwp) ::  jg    !<
    1333     INTEGER(iwp) ::  jh    !<
    1334     INTEGER(iwp) ::  jink  !<
    1335     INTEGER(iwp) ::  jump  !<
    1336     INTEGER(iwp) ::  k     !<
    1337     INTEGER(iwp) ::  kb    !<
    1338     INTEGER(iwp) ::  kc    !<
    1339     INTEGER(iwp) ::  kd    !<
    1340     INTEGER(iwp) ::  ke    !<
    1341     INTEGER(iwp) ::  kf    !<
    1342     INTEGER(iwp) ::  kstop !<
    1343     INTEGER(iwp) ::  l     !<
    1344     INTEGER(iwp) ::  m     !<
    1345 
    1346     REAL(wp) ::  a10(nfft) !<
    1347     REAL(wp) ::  a11(nfft) !<
    1348     REAL(wp) ::  a20(nfft) !<
    1349     REAL(wp) ::  a21(nfft) !<
    1350     REAL(wp) ::  b10(nfft) !<
    1351     REAL(wp) ::  b11(nfft) !<
    1352     REAL(wp) ::  b20(nfft) !<
    1353     REAL(wp) ::  b21(nfft) !<
     1314    INTEGER(iwp) ::  i      !<
     1315    INTEGER(iwp) ::  ia     !<
     1316    INTEGER(iwp) ::  ib     !<
     1317    INTEGER(iwp) ::  ibase  !<
     1318    INTEGER(iwp) ::  ic     !<
     1319    INTEGER(iwp) ::  id     !<
     1320    INTEGER(iwp) ::  ie     !<
     1321    INTEGER(iwp) ::  if     !<
     1322    INTEGER(iwp) ::  igo    !<
     1323    INTEGER(iwp) ::  iink   !<
     1324    INTEGER(iwp) ::  ijk    !<
     1325    INTEGER(iwp) ::  j      !<
     1326    INTEGER(iwp) ::  ja     !<
     1327    INTEGER(iwp) ::  jb     !<
     1328    INTEGER(iwp) ::  jbase  !<
     1329    INTEGER(iwp) ::  jc     !<
     1330    INTEGER(iwp) ::  jd     !<
     1331    INTEGER(iwp) ::  je     !<
     1332    INTEGER(iwp) ::  jf     !<
     1333    INTEGER(iwp) ::  jg     !<
     1334    INTEGER(iwp) ::  jh     !<
     1335    INTEGER(iwp) ::  jink   !<
     1336    INTEGER(iwp) ::  jump   !<
     1337    INTEGER(iwp) ::  k      !<
     1338    INTEGER(iwp) ::  kb     !<
     1339    INTEGER(iwp) ::  kc     !<
     1340    INTEGER(iwp) ::  kd     !<
     1341    INTEGER(iwp) ::  ke     !<
     1342    INTEGER(iwp) ::  kf     !<
     1343    INTEGER(iwp) ::  kstop  !<
     1344    INTEGER(iwp) ::  l      !<
     1345    INTEGER(iwp) ::  m      !<
     1346
     1347    REAL(wp) ::  a(*)      !<
     1348    REAL(wp) ::  b(*)      !<
     1349    REAL(wp) ::  c(*)      !<
     1350    REAL(wp) ::  d(*)      !<
     1351    REAL(wp) ::  trigs(*)  !<
     1352
     1353    REAL(wp) ::  c1      !<
     1354    REAL(wp) ::  c2      !<
     1355    REAL(wp) ::  c3      !<
     1356    REAL(wp) ::  c4      !<
     1357    REAL(wp) ::  c5      !<
     1358    REAL(wp) ::  qqrt5   !<
     1359    REAL(wp) ::  qrt5    !<
     1360    REAL(wp) ::  s1      !<
     1361    REAL(wp) ::  s2      !<
     1362    REAL(wp) ::  s3      !<
     1363    REAL(wp) ::  s4      !<
     1364    REAL(wp) ::  s5      !<
     1365    REAL(wp) ::  sin36   !<
     1366    REAL(wp) ::  sin45   !<
     1367    REAL(wp) ::  sin60   !<
     1368    REAL(wp) ::  sin72   !<
     1369    REAL(wp) ::  ssin36  !<
     1370    REAL(wp) ::  ssin45  !<
     1371    REAL(wp) ::  ssin60  !<
     1372    REAL(wp) ::  ssin72  !<
     1373
     1374    REAL(wp) ::  a10(nfft)  !<
     1375    REAL(wp) ::  a11(nfft)  !<
     1376    REAL(wp) ::  a20(nfft)  !<
     1377    REAL(wp) ::  a21(nfft)  !<
     1378    REAL(wp) ::  b10(nfft)  !<
     1379    REAL(wp) ::  b11(nfft)  !<
     1380    REAL(wp) ::  b20(nfft)  !<
     1381    REAL(wp) ::  b21(nfft)  !<
    13541382
    13551383
     
    13631391    iink = la * inc1
    13641392    jink = la * inc2
    1365     jump = (ifac-1) * jink
    1366     kstop = (n-ifac) / (2*ifac)
     1393    jump = ( ifac - 1 ) * jink
     1394    kstop = ( n - ifac ) / ( 2 * ifac )
    13671395
    13681396    IF ( lot > nfft )  THEN
     
    13861414
    13871415          ia = 1
    1388           ib = ia + (2*m-la) * inc1
     1416          ib = ia + ( 2 * m - la ) * inc1
    13891417          ja = 1
    13901418          jb = ja + jink
     
    14091437
    14101438             ia = ia + iink
    1411              iink = 2 * iink
     1439             iink  = 2 * iink
    14121440             ib = ib - iink
    14131441             ibase = 0
     
    14311459                         c(ja+j) = a(ia+i) + a(ib+i)
    14321460                         d(ja+j) = b(ia+i) - b(ib+i)
    1433                          c(jb+j) = c1*(a(ia+i)-a(ib+i)) - s1*(b(ia+i)+b(ib+i))
    1434                          d(jb+j) = s1*(a(ia+i)-a(ib+i)) + c1*(b(ia+i)+b(ib+i))
     1461                         c(jb+j) = c1 * ( a(ia+i) - a(ib+i) ) - s1 * ( b(ia+i) + b(ib+i) )
     1462                         d(jb+j) = s1 * ( a(ia+i) - a(ib+i) ) + c1 * ( b(ia+i) + b(ib+i) )
    14351463                         i = i + inc3
    14361464                         j = j + inc4
     
    14571485                DO  ijk = 1, lot
    14581486                   c(ja+j) = a(ia+i)
    1459                    c(jb+j) = -b(ia+i)
     1487                   c(jb+j) = - b(ia+i)
    14601488                   i = i + inc3
    14611489                   j = j + inc4
     
    14741502                !OCL NOVREC
    14751503                DO  ijk = 1, lot
    1476                    c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i))
    1477                    c(jb+j) = 2.0_wp*(a(ia+i)-a(ib+i))
     1504                   c(ja+j) = 2.0_wp * ( a(ia+i) + a(ib+i) )
     1505                   c(jb+j) = 2.0_wp * ( a(ia+i) - a(ib+i) )
    14781506                   i = i + inc3
    14791507                   j = j + inc4
     
    14901518
    14911519          ia = 1
    1492           ib = ia + (2*m-la) * inc1
     1520          ib = ia + ( 2 * m - la ) * inc1
    14931521          ic = ib
    14941522          ja = 1
     
    15061534                DO  ijk = 1, lot
    15071535                   c(ja+j) = a(ia+i) + a(ib+i)
    1508                    c(jb+j) = (a(ia+i)-0.5_wp*a(ib+i)) - (sin60*(b(ib+i)))
    1509                    c(jc+j) = (a(ia+i)-0.5_wp*a(ib+i)) + (sin60*(b(ib+i)))
     1536                   c(jb+j) = ( a(ia+i) - 0.5_wp * a(ib+i) ) - ( sin60 * ( b(ib+i) ) )
     1537                   c(jc+j) = ( a(ia+i) - 0.5_wp * a(ib+i) ) + ( sin60 * ( b(ib+i) ) )
    15101538                   i = i + inc3
    15111539                   j = j + inc4
     
    15381566                      !OCL NOVREC
    15391567                      DO  ijk = 1, lot
    1540                          c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
    1541                          d(ja+j) = b(ia+i) + (b(ib+i)-b(ic+i))
    1542                          c(jb+j) = c1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ b(ic+i)))) &
    1543                                    - s1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- a(ic+i))))
    1544                          d(jb+j) = s1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ b(ic+i)))) &
    1545                                    + c1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- a(ic+i))))
    1546                          c(jc+j) = c2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ b(ic+i)))) &
    1547                                    - s2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- a(ic+i))))
    1548                          d(jc+j) = s2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ b(ic+i)))) &
    1549                                    + c2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- a(ic+i))))
     1568                         c(ja+j) = a(ia+i) + ( a(ib+i) + a(ic+i) )
     1569                         d(ja+j) = b(ia+i) + ( b(ib+i) - b(ic+i) )
     1570                         c(jb+j) = c1 * ( ( a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) ) ) - ( sin60   &
     1571                                   * ( b(ib+i) + b(ic+i) ) ) ) - s1 * ( ( b(ia+i) - 0.5_wp         &
     1572                                   * ( b(ib+i) - b(ic+i) ) ) + ( sin60 * ( a(ib+i) - a(ic+i) ) ) )
     1573                         d(jb+j) = s1 * ( ( a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) ) ) - ( sin60   &
     1574                                   * ( b(ib+i) + b(ic+i) ) ) ) + c1 * ( ( b(ia+i) - 0.5_wp         &
     1575                                   * ( b(ib+i) - b(ic+i) ) ) + ( sin60 * ( a(ib+i) - a(ic+i) ) ) )
     1576                         c(jc+j) = c2 * ( ( a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) ) ) + ( sin60   &
     1577                                   * ( b(ib+i) + b(ic+i) ) ) ) - s2 * ( ( b(ia+i) - 0.5_wp         &
     1578                                   * ( b(ib+i) - b(ic+i) ) ) - ( sin60 * ( a(ib+i) - a(ic+i) ) ) )
     1579                         d(jc+j) = s2 * ( ( a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) ) ) + ( sin60   &
     1580                                   * ( b(ib+i) + b(ic+i) ) ) ) + c2 * ( ( b(ia+i) - 0.5_wp         &
     1581                                   * ( b(ib+i) - b(ic+i) ) ) - ( sin60 * ( a(ib+i) - a(ic+i) ) ) )
    15501582                         i = i + inc3
    15511583                         j = j + inc4
     
    15731605                DO  ijk = 1, lot
    15741606                   c(ja+j) = a(ia+i) + a(ib+i)
    1575                    c(jb+j) = (0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
    1576                    c(jc+j) = -(0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
     1607                   c(jb+j) = ( 0.5_wp * a(ia+i) - a(ib+i) ) - ( sin60 * b(ia+i) )
     1608                   c(jc+j) = - ( 0.5_wp * a(ia+i) - a(ib+i) ) - ( sin60 * b(ia+i) )
    15771609                   i = i + inc3
    15781610                   j = j + inc4
     
    15921624                !OCL NOVREC
    15931625                DO  ijk = 1, lot
    1594                    c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i))
    1595                    c(jb+j) = (2.0_wp*a(ia+i)-a(ib+i)) - (ssin60*b(ib+i))
    1596                    c(jc+j) = (2.0_wp*a(ia+i)-a(ib+i)) + (ssin60*b(ib+i))
     1626                   c(ja+j) = 2.0_wp * ( a(ia+i) + a(ib+i) )
     1627                   c(jb+j) = ( 2.0_wp * a(ia+i) - a(ib+i) ) - ( ssin60 * b(ib+i) )
     1628                   c(jc+j) = ( 2.0_wp * a(ia+i) - a(ib+i) ) + ( ssin60 * b(ib+i) )
    15971629                   i = i + inc3
    15981630                   j = j + inc4
     
    16091641
    16101642          ia = 1
    1611           ib = ia + (2*m-la) * inc1
    1612           ic = ib + 2*m*inc1
     1643          ib = ia + ( 2 * m - la ) * inc1
     1644          ic = ib + 2 * m * inc1
    16131645          id = ib
    16141646          ja = 1
     
    16261658                !OCL NOVREC
    16271659                DO  ijk = 1, lot
    1628                    c(ja+j) = (a(ia+i)+a(ic+i)) + a(ib+i)
    1629                    c(jb+j) = (a(ia+i)-a(ic+i)) - b(ib+i)
    1630                    c(jc+j) = (a(ia+i)+a(ic+i)) - a(ib+i)
    1631                    c(jd+j) = (a(ia+i)-a(ic+i)) + b(ib+i)
     1660                   c(ja+j) = ( a(ia+i) + a(ic+i) ) + a(ib+i)
     1661                   c(jb+j) = ( a(ia+i) - a(ic+i) ) - b(ib+i)
     1662                   c(jc+j) = ( a(ia+i) + a(ic+i) ) - a(ib+i)
     1663                   c(jd+j) = ( a(ia+i) - a(ic+i) ) + b(ib+i)
    16321664                   i = i + inc3
    16331665                   j = j + inc4
     
    16361668                jbase = jbase + inc2
    16371669             ENDDO
    1638              ia = ia + iink
    1639              iink = 2 * iink
    1640              ib = ib + iink
    1641              ic = ic - iink
    1642              id = id - iink
     1670             ia    = ia + iink
     1671             iink  = 2 * iink
     1672             ib    = ib + iink
     1673             ic    = ic - iink
     1674             id    = id - iink
    16431675             jbase = jbase + jump
    1644              jump = 2 * jump + jink
     1676             jump  = 2 * jump + jink
    16451677
    16461678             IF ( ib /= ic )  THEN
     
    16641696                      !OCL NOVREC
    16651697                      DO  ijk = 1, lot
    1666                          c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))
    1667                          d(ja+j) = (b(ia+i)-b(ic+i)) + (b(ib+i)-b(id+i))
    1668                          c(jc+j) = c2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) - s2*((b(ia+i)-b(ic+i))&
    1669                                    -(b(ib+i)-b(id+i)))
    1670                          d(jc+j) = s2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) + c2*((b(ia+i)-b(ic+i))&
    1671                                    -(b(ib+i)-b(id+i)))
    1672                          c(jb+j) = c1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) - s1*((b(ia+i)+b(ic+i))&
    1673                                    +(a(ib+i)-a(id+i)))
    1674                          d(jb+j) = s1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) + c1*((b(ia+i)+b(ic+i))&
    1675                                    +(a(ib+i)-a(id+i)))
    1676                          c(jd+j) = c3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) - s3*((b(ia+i)+b(ic+i))&
    1677                                    -(a(ib+i)-a(id+i)))
    1678                          d(jd+j) = s3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) + c3*((b(ia+i)+b(ic+i))&
    1679                                    -(a(ib+i)-a(id+i)))
     1698                         c(ja+j) = ( a(ia+i) + a(ic+i) ) + ( a(ib+i) + a(id+i) )
     1699                         d(ja+j) = ( b(ia+i) - b(ic+i) ) + ( b(ib+i) - b(id+i) )
     1700                         c(jc+j) = c2 * ( ( a(ia+i) + a(ic+i) ) - ( a(ib+i) + a(id+i) ) ) - s2     &
     1701                                   * ( ( b(ia+i) - b(ic+i) ) - ( b(ib+i) - b(id+i) ) )
     1702                         d(jc+j) = s2 * ( ( a(ia+i) + a(ic+i) ) - ( a(ib+i) + a(id+i) ) ) + c2     &
     1703                                   * ( ( b(ia+i) - b(ic+i) ) - ( b(ib+i) - b(id+i) ) )
     1704                         c(jb+j) = c1 * ( ( a(ia+i) - a(ic+i) ) - ( b(ib+i) + b(id+i) ) ) - s1     &
     1705                                   * ( ( b(ia+i) + b(ic+i) ) + ( a(ib+i) - a(id+i) ) )
     1706                         d(jb+j) = s1 * ( ( a(ia+i) - a(ic+i) ) - ( b (ib+i) + b(id+i) ) ) + c1    &
     1707                                   * ( ( b(ia+i) + b(ic+i) ) + ( a(ib+i) - a(id+i) ) )
     1708                         c(jd+j) = c3 * ( ( a(ia+i) - a(ic+i) ) + ( b(ib+i) + b(id+i) ) ) - s3     &
     1709                                   * ( ( b(ia+i) + b(ic+i) ) - ( a(ib+i) - a(id+i) ) )
     1710                         d(jd+j) = s3 * ( ( a(ia+i) - a(ic+i) ) + ( b(ib+i) + b(id+i) ) ) + c3     &
     1711                                   * ( ( b(ia+i) + b(ic+i) ) - ( a(ib+i) - a(id+i) ) )
    16801712                         i = i + inc3
    16811713                         j = j + inc4
     
    17051737                DO  ijk = 1, lot
    17061738                   c(ja+j) = a(ia+i) + a(ib+i)
    1707                    c(jb+j) = sin45*((a(ia+i)-a(ib+i))-(b(ia+i)+b(ib+i)))
     1739                   c(jb+j) = sin45 * ( ( a(ia+i) - a(ib+i) ) - ( b(ia+i) + b(ib+i) ) )
    17081740                   c(jc+j) = b(ib+i) - b(ia+i)
    1709                    c(jd+j) = -sin45*((a(ia+i)-a(ib+i))+(b(ia+i)+b(ib+i)))
     1741                   c(jd+j) = - sin45 * ( ( a(ia+i) - a(ib+i) ) + ( b(ia+i) + b(ib+i) ) )
    17101742                   i = i + inc3
    17111743                   j = j + inc4
     
    17241756                !OCL NOVREC
    17251757                DO  ijk = 1, lot
    1726                    c(ja+j) = 2.0_wp*((a(ia+i)+a(ic+i))+a(ib+i))
    1727                    c(jb+j) = 2.0_wp*((a(ia+i)-a(ic+i))-b(ib+i))
    1728                    c(jc+j) = 2.0_wp*((a(ia+i)+a(ic+i))-a(ib+i))
    1729                    c(jd+j) = 2.0_wp*((a(ia+i)-a(ic+i))+b(ib+i))
     1758                   c(ja+j) = 2.0_wp * ( ( a(ia+i) + a(ic+i) ) + a(ib+i) )
     1759                   c(jb+j) = 2.0_wp * ( ( a(ia+i) - a(ic+i) ) - b(ib+i) )
     1760                   c(jc+j) = 2.0_wp * ( ( a(ia+i) + a(ic+i) ) - a(ib+i) )
     1761                   c(jd+j) = 2.0_wp * ( ( a(ia+i) - a(ic+i) ) + b(ib+i) )
    17301762                   i = i + inc3
    17311763                   j = j + inc4
     
    17421774
    17431775          ia = 1
    1744           ib = ia + (2*m-la) * inc1
    1745           ic = ib + 2*m*inc1
     1776          ib = ia + ( 2 * m - la ) * inc1
     1777          ic = ib + 2 * m * inc1
    17461778          id = ic
    17471779          ie = ib
     
    17611793                !OCL NOVREC
    17621794                DO  ijk = 1, lot
    1763                    c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
    1764                    c(jb+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) -        &
    1765                              (sin72*b(ib+i)+sin36*b(ic+i))
    1766                    c(jc+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) -        &
    1767                              (sin36*b(ib+i)-sin72*b(ic+i))
    1768                    c(jd+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) +        &
    1769                              (sin36*b(ib+i)-sin72*b(ic+i))
    1770                    c(je+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) +        &
    1771                              (sin72*b(ib+i)+sin36*b(ic+i))
     1795                   c(ja+j) = a(ia+i) + ( a(ib+i) + a(ic+i) )
     1796                   c(jb+j) = ( ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) + qrt5 *              &
     1797                             ( a(ib+i) - a(ic+i) ) ) - ( sin72 * b(ib+i) + sin36 * b(ic+i) )
     1798                   c(jc+j) = ( ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) - qrt5 *              &
     1799                             ( a(ib+i) - a(ic+i) ) ) - ( sin36 * b(ib+i) - sin72 * b(ic+i) )
     1800                   c(jd+j) = ( ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) - qrt5 *              &
     1801                             ( a(ib+i) - a(ic+i) ) ) + ( sin36 * b(ib+i) - sin72 * b(ic+i) )
     1802                   c(je+j) = ( ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) + qrt5 *              &
     1803                             ( a(ib+i) - a(ic+i) ) ) + ( sin72 * b(ib+i) + sin36 * b(ic+i) )
    17721804                   i = i + inc3
    17731805                   j = j + inc4
     
    17761808                jbase = jbase + inc2
    17771809             ENDDO
    1778              ia = ia + iink
    1779              iink = 2 * iink
    1780              ib = ib + iink
    1781              ic = ic + iink
    1782              id = id - iink
    1783              ie = ie - iink
     1810             ia    = ia + iink
     1811             iink  = 2 * iink
     1812             ib    = ib + iink
     1813             ic    = ic + iink
     1814             id    = id - iink
     1815             ie    = ie - iink
    17841816             jbase = jbase + jump
    1785              jump = 2 * jump + jink
     1817             jump  = 2 * jump + jink
    17861818
    17871819             IF ( ib /= id )  THEN
     
    18081840                      !OCL NOVREC
    18091841                      DO  ijk = 1, lot
    1810                          a10(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) +      &
    1811                                     qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
    1812                          a20(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) -      &
    1813                                     qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
    1814                          b10(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) +      &
    1815                                     qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
    1816                          b20(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) -      &
    1817                                     qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
    1818                          a11(ijk) = sin72*(b(ib+i)+b(ie+i)) + sin36*(b(ic+i)+b(id+i))
    1819                          a21(ijk) = sin36*(b(ib+i)+b(ie+i)) - sin72*(b(ic+i)+b(id+i))
    1820                          b11(ijk) = sin72*(a(ib+i)-a(ie+i)) + sin36*(a(ic+i)-a(id+i))
    1821                          b21(ijk) = sin36*(a(ib+i)-a(ie+i)) - sin72*(a(ic+i)-a(id+i))
    1822                          c(ja+j) = a(ia+i) + ((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))
    1823                          d(ja+j) = b(ia+i) + ((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))
    1824                          c(jb+j) = c1*(a10(ijk)-a11(ijk)) - s1*(b10(ijk)+b11(ijk))
    1825                          d(jb+j) = s1*(a10(ijk)-a11(ijk)) + c1*(b10(ijk)+b11(ijk))
    1826                          c(je+j) = c4*(a10(ijk)+a11(ijk)) - s4*(b10(ijk)-b11(ijk))
    1827                          d(je+j) = s4*(a10(ijk)+a11(ijk)) + c4*(b10(ijk)-b11(ijk))
    1828                          c(jc+j) = c2*(a20(ijk)-a21(ijk)) - s2*(b20(ijk)+b21(ijk))
    1829                          d(jc+j) = s2*(a20(ijk)-a21(ijk)) + c2*(b20(ijk)+b21(ijk))
    1830                          c(jd+j) = c3*(a20(ijk)+a21(ijk)) - s3*(b20(ijk)-b21(ijk))
    1831                          d(jd+j) = s3*(a20(ijk)+a21(ijk)) + c3*(b20(ijk)-b21(ijk))
    1832 
     1842                         a10(ijk) = ( a(ia+i) - 0.25_wp * ( ( a(ib+i) + a(ie+i) ) + ( a(ic+i)      &
     1843                                    + a(id+i) ) ) ) + qrt5 * ( ( a(ib+i) + a(ie+i) ) - ( a(ic+i)   &
     1844                                    + a(id+i) ) )
     1845                         a20(ijk) = ( a(ia+i) - 0.25_wp * ( ( a(ib+i) + a(ie+i) ) + ( a(ic+i)      &
     1846                                    + a(id+i) ) ) ) - qrt5 * ( ( a(ib+i) + a(ie+i) ) - ( a(ic+i)   &
     1847                                    + a(id+i) ) )
     1848                         b10(ijk) = ( b(ia+i) - 0.25_wp * ( ( b(ib+i) - b(ie+i) ) + ( b(ic+i)      &
     1849                                    - b(id+i) ) ) ) + qrt5 * ( ( b(ib+i) - b(ie+i) ) - ( b(ic+i)   &
     1850                                    - b(id+i) ) )
     1851                         b20(ijk) = ( b(ia+i) - 0.25_wp * ( ( b(ib+i) - b(ie+i) ) + ( b(ic+i)      &
     1852                                    - b(id+i) ) ) ) - qrt5 * ( ( b(ib+i) - b(ie+i) ) - ( b(ic+i)   &
     1853                                    - b(id+i) ) )
     1854                         a11(ijk) = sin72 * ( b(ib+i) + b(ie+i) ) + sin36 * ( b(ic+i) + b(id+i) )
     1855                         a21(ijk) = sin36 * ( b(ib+i) + b(ie+i) ) - sin72 * ( b(ic+i) + b(id+i) )
     1856                         b11(ijk) = sin72 * ( a(ib+i) - a(ie+i) ) + sin36 * ( a(ic+i) - a(id+i) )
     1857                         b21(ijk) = sin36 * ( a(ib+i) - a(ie+i) ) - sin72 * ( a(ic+i) - a(id+i) )
     1858                         c(ja+j)  = a(ia+i) + ( ( a(ib+i) + a(ie+i) ) + ( a(ic+i) + a(id+i) ) )
     1859                         d(ja+j)  = b(ia+i) + ( ( b(ib+i) - b(ie+i) ) + ( b(ic+i) - b(id+i) ) )
     1860                         c(jb+j)  = c1 * ( a10(ijk) - a11(ijk) ) - s1 * ( b10(ijk) + b11(ijk) )
     1861                         d(jb+j)  = s1 * ( a10(ijk) - a11(ijk) ) + c1 * ( b10(ijk) + b11(ijk) )
     1862                         c(je+j)  = c4 * ( a10(ijk) + a11(ijk) ) - s4 * ( b10(ijk) - b11(ijk) )
     1863                         d(je+j)  = s4 * ( a10(ijk) + a11(ijk) ) + c4 * ( b10(ijk) - b11(ijk) )
     1864                         c(jc+j)  = c2 * ( a20(ijk) - a21(ijk) ) - s2 * ( b20(ijk) + b21(ijk) )
     1865                         d(jc+j)  = s2 * ( a20(ijk) - a21(ijk) ) + c2 * ( b20(ijk) + b21(ijk) )
     1866                         c(jd+j)  = c3 * ( a20(ijk) + a21(ijk) ) - s3 * ( b20(ijk) - b21(ijk) )
     1867                         d(jd+j)  = s3 * ( a20(ijk) + a21(ijk) ) + c3 * ( b20(ijk) - b21(ijk) )
    18331868                         i = i + inc3
    18341869                         j = j + inc4
     
    18571892                !OCL NOVREC
    18581893                DO  ijk = 1, lot
    1859                    c(ja+j) = (a(ia+i)+a(ib+i)) + a(ic+i)
    1860                    c(jb+j) = (qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) -        &
    1861                              (sin36*b(ia+i)+sin72*b(ib+i))
    1862                    c(je+j) = -(qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) -       &
    1863                              (sin36*b(ia+i)+sin72*b(ib+i))
    1864                    c(jc+j) = (qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) -        &
    1865                              (sin72*b(ia+i)-sin36*b(ib+i))
    1866                    c(jd+j) = -(qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) -       &
    1867                              (sin72*b(ia+i)-sin36*b(ib+i))
     1894                   c(ja+j) = ( a(ia+i) + a(ib+i) ) + a(ic+i)
     1895                   c(jb+j) = ( qrt5 * ( a(ia+i) - a(ib+i) ) + ( 0.25_wp * ( a(ia+i) + a(ib+i) )    &
     1896                             - a(ic+i) ) ) - ( sin36 * b(ia+i) + sin72 * b(ib+i) )
     1897                   c(je+j) = - ( qrt5 * ( a(ia+i) - a(ib+i) ) + ( 0.25_wp * ( a(ia+i) + a(ib+i) )  &
     1898                             - a(ic+i) ) ) - ( sin36 * b(ia+i) + sin72 * b(ib+i) )
     1899                   c(jc+j) = ( qrt5 * ( a(ia+i) - a(ib+i) ) - ( 0.25_wp * ( a(ia+i) + a(ib+i) )    &
     1900                             - a(ic+i) ) ) - ( sin72 * b(ia+i) - sin36 * b(ib+i) )
     1901                   c(jd+j) = - ( qrt5 * ( a(ia+i) - a(ib+i) ) - ( 0.25_wp * ( a(ia+i) + a(ib+i) )  &
     1902                             - a(ic+i) ) ) - ( sin72 * b(ia+i) - sin36 * b(ib+i) )
    18681903                   i = i + inc3
    18691904                   j = j + inc4
     
    18851920                !OCL NOVREC
    18861921                DO  ijk = 1, lot
    1887                    c(ja+j) = 2.0_wp*(a(ia+i)+(a(ib+i)+a(ic+i)))
    1888                    c(jb+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+i)))  &
    1889                              - (ssin72*b(ib+i)+ssin36*b(ic+i))
    1890                    c(jc+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+i)))  &
    1891                              - (ssin36*b(ib+i)-ssin72*b(ic+i))
    1892                    c(jd+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+i)))  &
    1893                              + (ssin36*b(ib+i)-ssin72*b(ic+i))
    1894                    c(je+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+i)))  &
    1895                              + (ssin72*b(ib+i)+ssin36*b(ic+i))
     1922                   c(ja+j) = 2.0_wp * ( a(ia+i) + ( a(ib+i) + a(ic+i) ) )
     1923                   c(jb+j) = ( 2.0_wp * ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) + qqrt5      &
     1924                             * ( a(ib+i) - a(ic+i) ) ) - ( ssin72 * b(ib+i) + ssin36 * b(ic+i) )
     1925                   c(jc+j) = ( 2.0_wp * ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) - qqrt5      &
     1926                             * ( a(ib+i) - a(ic+i) ) ) - ( ssin36 * b(ib+i) - ssin72 * b(ic+i) )
     1927                   c(jd+j) = ( 2.0_wp * ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) - qqrt5      &
     1928                             * ( a(ib+i) - a(ic+i) ) ) + ( ssin36 * b(ib+i) - ssin72 * b(ic+i) )
     1929                   c(je+j) = ( 2.0_wp * ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) + qqrt5      &
     1930                             * ( a(ib+i) - a(ic+i) ) ) + (ssin72 * b(ib+i) + ssin36 * b(ic+i) )
    18961931                   i = i + inc3
    18971932                   j = j + inc4
     
    19081943
    19091944          ia = 1
    1910           ib = ia + (2*m-la) * inc1
    1911           ic = ib + 2*m*inc1
    1912           id = ic + 2*m*inc1
     1945          ib = ia + ( 2 * m - la ) * inc1
     1946          ic = ib + 2 * m * inc1
     1947          id = ic + 2 * m * inc1
    19131948          ie = ic
    19141949          if = ib
     
    19291964                !OCL NOVREC
    19301965                DO  ijk = 1, lot
    1931                    c(ja+j) = (a(ia+i)+a(id+i)) + (a(ib+i)+a(ic+i))
    1932                    c(jd+j) = (a(ia+i)-a(id+i)) - (a(ib+i)-a(ic+i))
    1933                    c(jb+j) = ((a(ia+i)-a(id+i))+0.5_wp*(a(ib+i)-a(ic+i))) - (sin60*(b(ib+i)+b(ic+i)))
    1934                    c(jf+j) = ((a(ia+i)-a(id+i))+0.5_wp*(a(ib+i)-a(ic+i))) + (sin60*(b(ib+i)+b(ic+i)))
    1935                    c(jc+j) = ((a(ia+i)+a(id+i))-0.5_wp*(a(ib+i)+a(ic+i))) - (sin60*(b(ib+i)-b(ic+i)))
    1936                    c(je+j) = ((a(ia+i)+a(id+i))-0.5_wp*(a(ib+i)+a(ic+i))) + (sin60*(b(ib+i)-b(ic+i)))
     1966                   c(ja+j) = ( a (ia+i) + a(id+i) ) +  ( a(ib+i) + a(ic+i) )
     1967                   c(jd+j) = ( a (ia+i) - a(id+i) ) -  ( a(ib+i) - a(ic+i) )
     1968                   c(jb+j) = ( ( a(ia+i) - a(id+i) ) + 0.5_wp * ( a(ib+i) - a(ic+i) ) )            &
     1969                             - ( sin60 * ( b(ib+i) + b(ic+i) ) )
     1970                   c(jf+j) = ( ( a(ia+i) - a(id+i) ) + 0.5_wp * ( a(ib+i) - a(ic+i) ) )            &
     1971                             + ( sin60 * b(ib+i) + b(ic+i) )
     1972                   c(jc+j) = ( ( a(ia+i) + a(id+i) ) - 0.5_wp * ( a(ib+i) + a(ic+i) ) )            &
     1973                             - ( sin60 * b(ib+i) - b(ic+i) )
     1974                   c(je+j) = ( ( a(ia+i) + a(id+i) ) - 0.5_wp * ( a(ib+i) + a(ic+i) ) )            &
     1975                             + ( sin60 * b(ib+i) - b(ic+i) )
    19371976                   i = i + inc3
    19381977                   j = j + inc4
     
    19411980                jbase = jbase + inc2
    19421981             ENDDO
    1943              ia = ia + iink
    1944              iink = 2 * iink
    1945              ib = ib + iink
    1946              ic = ic + iink
    1947              id = id - iink
    1948              ie = ie - iink
    1949              if = if - iink
     1982             ia    = ia + iink
     1983             iink  = 2 * iink
     1984             ib    = ib + iink
     1985             ic    = ic + iink
     1986             id    = id - iink
     1987             ie    = ie - iink
     1988             if    = if - iink
    19501989             jbase = jbase + jump
    1951              jump = 2 * jump + jink
     1990             jump  = 2 * jump + jink
    19521991
    19531992             IF ( ic /= id )  THEN
     
    19772016                      !OCL NOVREC
    19782017                      DO  ijk = 1, lot
    1979                          a11(ijk) = (a(ie+i)+a(ib+i)) + (a(ic+i)+a(if+i))
    1980                          a20(ijk) = (a(ia+i)+a(id+i)) - 0.5_wp*a11(ijk)
    1981                          a21(ijk) = sin60*((a(ie+i)+a(ib+i))-(a(ic+i)+a(if+i)))
    1982                          b11(ijk) = (b(ib+i)-b(ie+i)) + (b(ic+i)-b(if+i))
    1983                          b20(ijk) = (b(ia+i)-b(id+i)) - 0.5_wp*b11(ijk)
    1984                          b21(ijk) = sin60*((b(ib+i)-b(ie+i))-(b(ic+i)-b(if+i)))
    1985 
    1986                          c(ja+j) = (a(ia+i)+a(id+i)) + a11(ijk)
    1987                          d(ja+j) = (b(ia+i)-b(id+i)) + b11(ijk)
    1988                          c(jc+j) = c2*(a20(ijk)-b21(ijk)) - s2*(b20(ijk)+a21(ijk))
    1989                          d(jc+j) = s2*(a20(ijk)-b21(ijk)) + c2*(b20(ijk)+a21(ijk))
    1990                          c(je+j) = c4*(a20(ijk)+b21(ijk)) - s4*(b20(ijk)-a21(ijk))
    1991                          d(je+j) = s4*(a20(ijk)+b21(ijk)) + c4*(b20(ijk)-a21(ijk))
    1992 
    1993                          a11(ijk) = (a(ie+i)-a(ib+i)) + (a(ic+i)-a(if+i))
    1994                          b11(ijk) = (b(ie+i)+b(ib+i)) - (b(ic+i)+b(if+i))
    1995                          a20(ijk) = (a(ia+i)-a(id+i)) - 0.5_wp*a11(ijk)
    1996                          a21(ijk) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
    1997                          b20(ijk) = (b(ia+i)+b(id+i)) + 0.5_wp*b11(ijk)
    1998                          b21(ijk) = sin60*((b(ie+i)+b(ib+i))+(b(ic+i)+b(if+i)))
    1999 
    2000                          c(jd+j) = c3*((a(ia+i)-a(id+i))+a11(ijk)) - s3*((b(ia+i)+b(id+i))-b11(ijk))
    2001                          d(jd+j) = s3*((a(ia+i)-a(id+i))+a11(ijk)) + c3*((b(ia+i)+b(id+i))-b11(ijk))
    2002                          c(jb+j) = c1*(a20(ijk)-b21(ijk)) - s1*(b20(ijk)-a21(ijk))
    2003                          d(jb+j) = s1*(a20(ijk)-b21(ijk)) + c1*(b20(ijk)-a21(ijk))
    2004                          c(jf+j) = c5*(a20(ijk)+b21(ijk)) - s5*(b20(ijk)+a21(ijk))
    2005                          d(jf+j) = s5*(a20(ijk)+b21(ijk)) + c5*(b20(ijk)+a21(ijk))
     2018                         a11(ijk) = ( a(ie+i) + a(ib+i) ) + ( a(ic+i) + a(if+i) )
     2019                         a20(ijk) = ( a(ia+i) + a(id+i) ) - 0.5_wp * a11(ijk)
     2020                         a21(ijk) = sin60 * ( ( a(ie+i) + a(ib+i) ) - ( a(ic+i) + a(if+i) ) )
     2021                         b11(ijk) = ( b(ib+i) - b(ie+i) ) + ( b(ic+i) - b(if+i) )
     2022                         b20(ijk) = ( b(ia+i) - b(id+i) ) - 0.5_wp * b11(ijk)
     2023                         b21(ijk) = sin60 * ( ( b(ib+i) - b(ie+i) ) - ( b(ic+i) - b(if+i) ) )
     2024
     2025                         c(ja+j) = ( a(ia+i) + a(id+i) ) + a11(ijk)
     2026                         d(ja+j) = ( b(ia+i) - b(id+i) ) + b11(ijk)
     2027                         c(jc+j) = c2 * ( a20(ijk) - b21(ijk) ) - s2 * ( b20(ijk) + a21(ijk) )
     2028                         d(jc+j) = s2 * ( a20(ijk) - b21(ijk) ) + c2 * ( b20(ijk) + a21(ijk) )
     2029                         c(je+j) = c4 * ( a20(ijk) + b21(ijk) ) - s4 * ( b20(ijk) - a21(ijk) )
     2030                         d(je+j) = s4 * ( a20(ijk) + b21(ijk) ) + c4 * ( b20(ijk) - a21(ijk) )
     2031
     2032                         a11(ijk) = ( a(ie+i) - a(ib+i) ) + ( a(ic+i) - a(if+i) )
     2033                         b11(ijk) = ( b(ie+i) + b(ib+i) ) - ( b(ic+i) + b(if+i) )
     2034                         a20(ijk) = ( a(ia+i) - a(id+i) ) - 0.5_wp * a11(ijk)
     2035                         a21(ijk) = sin60 * ( ( a(ie+i) - a(ib+i) ) - ( a(ic+i) - a(if+i) ) )
     2036                         b20(ijk) = ( b(ia+i) + b(id+i) ) + 0.5_wp * b11(ijk)
     2037                         b21(ijk) = sin60 * ( ( b(ie+i) + b(ib+i) ) + ( b(ic+i) + b(if+i) ) )
     2038
     2039                         c(jd+j) = c3 * ( (a(ia+i) - a(id+i) ) + a11(ijk) ) - s3                   &
     2040                                      * ( (b(ia+i) + b(id+i) ) - b11(ijk) )
     2041                         d(jd+j) = s3 * ( (a(ia+i) - a(id+i) ) + a11(ijk) ) + c3                   &
     2042                                      * ( (b(ia+i) + b(id+i) ) - b11(ijk) )
     2043                         c(jb+j) = c1 * ( a20(ijk) - b21(ijk) ) - s1 * ( b20(ijk) - a21(ijk) )
     2044                         d(jb+j) = s1 * ( a20(ijk) - b21(ijk) ) + c1 * ( b20(ijk) - a21(ijk) )
     2045                         c(jf+j) = c5 * ( a20(ijk) + b21(ijk) ) - s5 * ( b20(ijk) + a21(ijk) )
     2046                         d(jf+j) = s5 * ( a20(ijk) + b21(ijk) ) + c5 * ( b20(ijk) + a21(ijk) )
    20062047
    20072048                         i = i + inc3
     
    20322073                !OCL NOVREC
    20332074                DO  ijk = 1, lot
    2034                    c(ja+j) = a(ib+i) + (a(ia+i)+a(ic+i))
    2035                    c(jd+j) = b(ib+i) - (b(ia+i)+b(ic+i))
    2036                    c(jb+j) = (sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i))
    2037                    c(jf+j) = -(sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i))
    2038                    c(jc+j) = sin60*(b(ic+i)-b(ia+i)) + (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i))
    2039                    c(je+j) = sin60*(b(ic+i)-b(ia+i)) - (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i))
     2075                   c(ja+j) = a(ib+i) + ( a(ia+i) + a(ic+i) )
     2076                   c(jd+j) = b(ib+i) - ( b(ia+i) + b(ic+i) )
     2077                   c(jb+j) = ( sin60 * ( a(ia+i) - a(ic+i) ) ) - ( 0.5_wp * ( b(ia+i) + b(ic+i) )  &
     2078                             + b(ib+i) )
     2079                   c(jf+j) = - ( sin60 * ( a(ia+i) - a(ic+i) ) ) - ( 0.5_wp * ( b(ia+i)            &
     2080                             + b(ic+i) ) + b(ib+i) )
     2081                   c(jc+j) = sin60 * ( b(ic+i) - b(ia+i) ) + ( 0.5_wp * ( a(ia+i) + a(ic+i) )      &
     2082                             - a(ib+i) )
     2083                   c(je+j) = sin60 * ( b(ic+i) - b(ia+i) ) - ( 0.5_wp * ( a(ia+i) + a(ic+i) )      &
     2084                             - a(ib+i) )
    20402085                   i = i + inc3
    20412086                   j = j + inc4
     
    20552100                !OCL NOVREC
    20562101                DO  ijk = 1, lot
    2057                    c(ja+j) = (2.0_wp*(a(ia+i)+a(id+i))) + (2.0_wp*(a(ib+i)+a(ic+i)))
    2058                    c(jd+j) = (2.0_wp*(a(ia+i)-a(id+i))) - (2.0_wp*(a(ib+i)-a(ic+i)))
    2059                    c(jb+j) = (2.0_wp*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) - (ssin60*(b(ib+i)+b(ic+i)))
    2060                    c(jf+j) = (2.0_wp*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) + (ssin60*(b(ib+i)+b(ic+i)))
    2061                    c(jc+j) = (2.0_wp*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) - (ssin60*(b(ib+i)-b(ic+i)))
    2062                    c(je+j) = (2.0_wp*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) + (ssin60*(b(ib+i)-b(ic+i)))
     2102                   c(ja+j) = ( 2.0_wp * ( a(ia+i) + a(id+i) ) ) + ( 2.0_wp * ( a(ib+i) + a(ic+i) ) )
     2103                   c(jd+j) = ( 2.0_wp * ( a(ia+i) - a(id+i) ) ) - ( 2.0_wp * ( a(ib+i) - a(ic+i) ) )
     2104                   c(jb+j) = ( 2.0_wp * ( a(ia+i) - a(id+i) ) + ( a(ib+i) - a(ic+i) ) )            &
     2105                             - ( ssin60 * ( b(ib+i) + b(ic+i) ) )
     2106                   c(jf+j) = ( 2.0_wp * ( a(ia+i) - a(id+i) ) + ( a(ib+i) - a(ic+i) ) )            &
     2107                             + ( ssin60 * ( b(ib+i) + b(ic+i) ) )
     2108                   c(jc+j) = ( 2.0_wp * ( a(ia+i) + a(id+i) ) - ( a(ib+i) + a(ic+i) ) )            &
     2109                             - ( ssin60 * ( b(ib+i) - b(ic+i) ) )
     2110                   c(je+j) = ( 2.0_wp * ( a(ia+i) + a(id+i) ) - ( a(ib+i) + a(ic+i) ) )            &
     2111                             + ( ssin60 * ( b(ib+i) - b(ic+i) ) )
    20632112                   i = i + inc3
    20642113                   j = j + inc4
     
    20802129
    20812130          ia = 1
    2082           ib = ia + la*inc1
    2083           ic = ib + 2*la*inc1
    2084           id = ic + 2*la*inc1
    2085           ie = id + 2*la*inc1
     2131          ib = ia + la * inc1
     2132          ic = ib + 2 * la * inc1
     2133          id = ic + 2 * la * inc1
     2134          ie = id + 2 * la * inc1
    20862135          ja = 1
    20872136          jb = ja + jink
     
    21012150             !OCL NOVREC
    21022151             DO  ijk = 1, lot
    2103                 c(ja+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))+(a(ib+i)+a(id+i)))
    2104                 c(je+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))-(a(ib+i)+a(id+i)))
    2105                 c(jc+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))-(b(ib+i)-b(id+i)))
    2106                 c(jg+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))+(b(ib+i)-b(id+i)))
    2107                 c(jb+j) = 2.0_wp*((a(ia+i)-a(ie+i))-b(ic+i)) + ssin45*((a(ib+i)-a(id+i))-(b(ib+i)+b(id+i)))
    2108                 c(jf+j) = 2.0_wp*((a(ia+i)-a(ie+i))-b(ic+i)) - ssin45*((a(ib+i)-a(id+i))-(b(ib+i)+b(id+i)))
    2109                 c(jd+j) = 2.0_wp*((a(ia+i)-a(ie+i))+b(ic+i)) - ssin45*((a(ib+i)-a(id+i))+(b(ib+i)+b(id+i)))
    2110                 c(jh+j) = 2.0_wp*((a(ia+i)-a(ie+i))+b(ic+i)) + ssin45*((a(ib+i)-a(id+i))+(b(ib+i)+b(id+i)))
     2152                c(ja+j) = 2.0_wp * ( ( ( a(ia+i) + a(ie+i) ) + a(ic+i) ) + ( a(ib+i) + a(id+i) ) )
     2153                c(je+j) = 2.0_wp * ( ( ( a(ia+i) + a(ie+i) ) + a(ic+i) ) - ( a(ib+i) + a(id+i) ) )
     2154                c(jc+j) = 2.0_wp * ( ( ( a(ia+i) + a(ie+i) ) - a(ic+i) ) - ( b(ib+i) - b(id+i) ) )
     2155                c(jg+j) = 2.0_wp * ( ( ( a(ia+i) + a(ie+i) ) - a(ic+i) ) + ( b(ib+i) - b(id+i) ) )
     2156                c(jb+j) = 2.0_wp * ( ( a(ia+i) - a(ie+i) ) - b(ic+i) ) + ssin45 * ( ( a(ib+i)      &
     2157                          - a(id+i) ) - ( b(ib+i) + b(id+i) ) )
     2158                c(jf+j) = 2.0_wp * ( ( a(ia+i) - a(ie+i) ) - b(ic+i) ) - ssin45 * ( ( a(ib+i)      &
     2159                          - a(id+i) ) - ( b(ib+i) + b(id+i) ) )
     2160                c(jd+j) = 2.0_wp * ( ( a(ia+i) - a(ie+i) ) + b(ic+i) ) - ssin45 * ( ( a(ib+i)      &
     2161                          - a(id+i) ) + ( b(ib+i) + b(id+i) ) )
     2162                c(jh+j) = 2.0_wp * ( ( a(ia+i) - a(ie+i) ) + b(ic+i) ) + ssin45 * ( ( a(ib+i)      &
     2163                          - a(id+i) ) + ( b(ib+i) + b(id+i) ) )
    21112164                i = i + inc3
    21122165                j = j + inc4
     
    21202173 END SUBROUTINE rpassm
    21212174
    2122 !------------------------------------------------------------------------------!
     2175!--------------------------------------------------------------------------------------------------!
    21232176! Description:
    21242177! ------------
    21252178!> Computes factors of n & trigonometric functins required by fft99 & fft991cy
    21262179!> Method: Dimension trigs(n),ifax(1),jfax(10),lfax(7)
    2127 !> subroutine 'set99' - computes factors of n & trigonometric
    2128 !> functins required by fft99 & fft991cy
    2129 !------------------------------------------------------------------------------!
     2180!> Subroutine 'set99' - computes factors of n & trigonometric functins required by fft99 & fft991cy
     2181!--------------------------------------------------------------------------------------------------!
    21302182 SUBROUTINE set99( trigs, ifax, n )
    21312183
     
    21382190    IMPLICIT NONE
    21392191
    2140     INTEGER(iwp) ::  n !<
    2141 
    2142     INTEGER(iwp) ::  ifax(*)  !<
    2143     REAL(wp)     ::  trigs(*) !<
    2144 
    2145 
    2146     REAL(wp) ::  angle    !<
    2147     REAL(wp) ::  del      !<
    2148     INTEGER(iwp) ::  i    !<
    2149     INTEGER(iwp) ::  ifac !<
    2150     INTEGER(iwp) ::  ixxx !<
    2151     INTEGER(iwp) ::  k    !<
    2152     INTEGER(iwp) ::  l    !<
    2153     INTEGER(iwp) ::  nfax !<
    2154     INTEGER(iwp) ::  nhl  !<
    2155     INTEGER(iwp) ::  nil  !<
    2156     INTEGER(iwp) ::  nu   !<
    2157 
    2158     INTEGER(iwp) ::  jfax(10) !<
    2159     INTEGER(iwp) ::  lfax(7)  !<
     2192    INTEGER(iwp) ::  ifax(*)   !<
     2193
     2194    INTEGER(iwp) ::  i         !<
     2195    INTEGER(iwp) ::  ifac      !<
     2196    INTEGER(iwp) ::  ixxx      !<
     2197    INTEGER(iwp) ::  jfax(10)  !<
     2198    INTEGER(iwp) ::  k         !<
     2199    INTEGER(iwp) ::  l         !<
     2200    INTEGER(iwp) ::  lfax(7)   !<
     2201    INTEGER(iwp) ::  n         !<
     2202    INTEGER(iwp) ::  nfax      !<
     2203    INTEGER(iwp) ::  nhl       !<
     2204    INTEGER(iwp) ::  nil       !<
     2205    INTEGER(iwp) ::  nu        !<
     2206
     2207    REAL(wp) ::  angle         !<
     2208    REAL(wp) ::  del           !<
     2209
     2210    REAL(wp) ::  trigs(*)      !<
    21602211
    21612212    DATA  lfax / 6, 8, 5, 4, 3, 2, 1 /
     
    21642215    ixxx = 1
    21652216
    2166     del = 4.0_wp * ASIN( 1.0_wp ) / REAL( n, KIND=wp )
     2217    del = 4.0_wp * ASIN( 1.0_wp ) / REAL( n, KIND = wp )
    21672218    nil = 0
    2168     nhl = (n/2) - 1
     2219    nhl = ( n / 2 ) - 1
    21692220    DO  k = nil, nhl
    2170        angle = REAL( k, KIND=wp ) * del
     2221       angle = REAL( k, KIND = wp ) * del
    21712222       trigs(2*k+1) = COS( angle )
    21722223       trigs(2*k+2) = SIN( angle )
     
    21982249       ifac = lfax(l)
    21992250       IF (ifac < 2 )  THEN
    2200           message_string = 'number of gridpoints along x or/and y ' //                             &
    2201                            'contain illegal  factors' //                                           &
     2251          message_string = 'number of gridpoints along x or/and y contain illegal  factors' //     &
    22022252                           '&only factors 2, 3, 5 are allowed'
    22032253          CALL message( 'temperton_fft', 'PA0311', 1, 2, 0, 6, 0 )
     
    22242274    IMPLICIT NONE
    22252275
    2226     REAL(wp),DIMENSION(:,:)     ::  a    !<
    2227     REAL(wp),DIMENSION(:)       ::  trigs !<
    2228     REAL(wp),DIMENSION(:,:)     ::  work  !<
     2276
     2277    INTEGER(iwp) ::  i       !<
     2278    INTEGER(iwp) ::  ia      !<
     2279    INTEGER(iwp) ::  ibase   !<
     2280    INTEGER(iwp) ::  ierr    !<
     2281    INTEGER(iwp) ::  ifac    !<
     2282    INTEGER(iwp) ::  igo     !<
     2283    INTEGER(iwp) ::  ii      !<
     2284    INTEGER(iwp) ::  inc     !<
     2285    INTEGER(iwp) ::  isign   !<
     2286    INTEGER(iwp) ::  istart  !<
     2287    INTEGER(iwp) ::  ix      !<
     2288    INTEGER(iwp) ::  iz      !<
     2289    INTEGER(iwp) ::  j       !<
     2290    INTEGER(iwp) ::  jbase   !<
     2291    INTEGER(iwp) ::  jump    !<
     2292    INTEGER(iwp) ::  k       !<
     2293    INTEGER(iwp) ::  la      !<
     2294    INTEGER(iwp) ::  lot     !<
     2295    INTEGER(iwp) ::  n       !<
     2296    INTEGER(iwp) ::  nfax    !<
     2297    INTEGER(iwp) ::  nvex    !<
     2298    INTEGER(iwp) ::  nx      !<
     2299
    22292300    INTEGER(iwp),DIMENSION(:),INTENT(IN) ::  ifax  !<
    22302301
    2231     INTEGER(iwp) ::  inc   !<
    2232     INTEGER(iwp) ::  isign !<
    2233     INTEGER(iwp) ::  jump  !<
    2234     INTEGER(iwp) ::  lot   !<
    2235     INTEGER(iwp) ::  n     !<
    2236 
    2237 
    2238     INTEGER(iwp) ::  i      !<
    2239     INTEGER(iwp) ::  ia     !<
    2240     INTEGER(iwp) ::  ibase  !<
    2241     INTEGER(iwp) ::  ierr   !<
    2242     INTEGER(iwp) ::  ifac   !<
    2243     INTEGER(iwp) ::  igo    !<
    2244     INTEGER(iwp) ::  ii     !<
    2245     INTEGER(iwp) ::  istart !<
    2246     INTEGER(iwp) ::  ix     !<
    2247     INTEGER(iwp) ::  iz     !<
    2248     INTEGER(iwp) ::  j      !<
    2249     INTEGER(iwp) ::  jbase  !<
    2250     INTEGER(iwp) ::  k      !<
    2251     INTEGER(iwp) ::  la     !<
    2252     INTEGER(iwp) ::  nfax   !<
    2253     INTEGER(iwp) ::  nvex   !<
    2254     INTEGER(iwp) ::  nx     !<
     2302    REAL(wp),DIMENSION(:,:) ::  a      !<
     2303    REAL(wp),DIMENSION(:)   ::  trigs  !<
     2304    REAL(wp),DIMENSION(:,:) ::  work   !<
    22552305
    22562306
     
    22622312    nfax = ifax(1)
    22632313    nx   = n + 1
    2264     IF ( MOD(n,2) == 1 )  nx = n
     2314    IF ( MOD( n, 2 ) == 1 )  nx = n
    22652315    nvex = 1
    22662316
     
    22732323       i = istart
    22742324       a(:,i+inc) = 0.5_wp * a(:,i)
    2275        IF ( MOD(n,2) /= 1 )  THEN
     2325       IF ( MOD( n, 2 ) /= 1 )  THEN
    22762326          i = istart + n * inc
    22772327          a(:,i) = 0.5_wp * a(:,i)
     
    23132363!
    23142364!--    If necessary, copy results back to a
    2315        IF ( MOD(nfax,2) /= 0 )  THEN
     2365       IF ( MOD( nfax, 2 ) /= 0 )  THEN
    23162366          ibase = 1
    23172367          jbase = ia
     
    23262376!
    23272377!--    Fill in zeros at end
    2328        ix = istart + n*inc
     2378       ix = istart + n * inc
    23292379       a(:,ix) = 0.0_wp
    23302380       a(:,ix+inc) = 0.0_wp
     
    23702420!
    23712421!--    If necessary, copy results back to a
    2372        IF ( MOD(nfax,2) /= 0 )  THEN
     2422       IF ( MOD( nfax, 2 ) /= 0 )  THEN
    23732423          ibase = 1
    23742424          jbase = ia
     
    23872437       a(:,ix+inc) = 0.0_wp
    23882438
    2389        IF ( MOD(n,2) /= 1 )  THEN
    2390           iz = istart + (n+1) * inc
     2439       IF ( MOD( n, 2 ) /= 1 )  THEN
     2440          iz = istart + ( n + 1 ) * inc
    23912441          a(:,iz) = 0.0_wp
    23922442       ENDIF
     
    23962446 END SUBROUTINE fft991cy_vec
    23972447
    2398 !------------------------------------------------------------------------------!
     2448!--------------------------------------------------------------------------------------------------!
    23992449! Description:
    24002450! ------------
    2401 !> Performs one pass through data as part of
    2402 !> multiple real fft (fourier analysis) routine.
     2451!> Performs one pass through data as part of multiple real fft (fourier analysis) routine.
    24032452!>
    24042453!> Method:
     
    24222471!>         2 - ifac not catered for
    24232472!>         3 - ifac only catered for if la=n/ifac
    2424 !------------------------------------------------------------------------------!
     2473!--------------------------------------------------------------------------------------------------!
    24252474 SUBROUTINE qpassm_vec( a, b, c, d, trigs, inc1, inc2, n, ifac, la, ierr )
    24262475
     
    24292478    IMPLICIT NONE
    24302479
    2431     INTEGER(iwp),INTENT(IN)  ::  ifac !<
    2432     INTEGER(iwp),INTENT(IN)  ::  inc1 !<
    2433     INTEGER(iwp),INTENT(IN)  ::  inc2 !<
    2434     INTEGER(iwp),INTENT(IN)  ::  la   !<
    2435     INTEGER(iwp),INTENT(IN)  ::  n    !<
    2436     INTEGER(iwp),INTENT(OUT) ::  ierr !<
     2480    INTEGER(iwp),INTENT(IN)  ::  ifac  !<
     2481    INTEGER(iwp),INTENT(IN)  ::  inc1  !<
     2482    INTEGER(iwp),INTENT(IN)  ::  inc2  !<
     2483    INTEGER(iwp),INTENT(IN)  ::  la    !<
     2484    INTEGER(iwp),INTENT(IN)  ::  n     !<
     2485
     2486    INTEGER(iwp),INTENT(OUT) ::  ierr  !<
     2487
     2488    INTEGER(iwp) ::  i      !<
     2489    INTEGER(iwp) ::  ia     !<
     2490    INTEGER(iwp) ::  ib     !<
     2491    INTEGER(iwp) ::  ibase  !<
     2492    INTEGER(iwp) ::  ic     !<
     2493    INTEGER(iwp) ::  id     !<
     2494    INTEGER(iwp) ::  ie     !<
     2495    INTEGER(iwp) ::  if     !<
     2496    INTEGER(iwp) ::  ig     !<
     2497    INTEGER(iwp) ::  igo    !<
     2498    INTEGER(iwp) ::  ih     !<
     2499    INTEGER(iwp) ::  iink   !<
     2500    INTEGER(iwp) ::  ijump  !<
     2501    INTEGER(iwp) ::  ipl    !<  loop index parallel loop
     2502    INTEGER(iwp) ::  j      !<
     2503    INTEGER(iwp) ::  ja     !<
     2504    INTEGER(iwp) ::  jb     !<
     2505    INTEGER(iwp) ::  jbase  !<
     2506    INTEGER(iwp) ::  jc     !<
     2507    INTEGER(iwp) ::  jd     !<
     2508    INTEGER(iwp) ::  je     !<
     2509    INTEGER(iwp) ::  jf     !<
     2510    INTEGER(iwp) ::  jink   !<
     2511    INTEGER(iwp) ::  k      !<
     2512    INTEGER(iwp) ::  kb     !<
     2513    INTEGER(iwp) ::  kc     !<
     2514    INTEGER(iwp) ::  kd     !<
     2515    INTEGER(iwp) ::  ke     !<
     2516    INTEGER(iwp) ::  kf     !<
     2517    INTEGER(iwp) ::  kstop  !<
     2518    INTEGER(iwp) ::  l      !<
     2519    INTEGER(iwp) ::  m      !<
    24372520
    24382521!
    24392522!-- Arrays are dimensioned with n
    2440     REAL(wp),DIMENSION(:,:) ::  a     !<
    2441     REAL(wp),DIMENSION(:,:) ::  b     !<
    2442     REAL(wp),DIMENSION(:,:) ::  c     !<
    2443     REAL(wp),DIMENSION(:,:) ::  d     !<
    2444     REAL(wp),DIMENSION(:),INTENT(IN) ::  trigs !<
    2445 
    2446     REAL(wp) ::  a0     !<
    2447     REAL(wp) ::  a1     !<
    2448     REAL(wp) ::  a10    !<
    2449     REAL(wp) ::  a11    !<
    2450     REAL(wp) ::  a2     !<
    2451     REAL(wp) ::  a20    !<
    2452     REAL(wp) ::  a21    !<
    2453     REAL(wp) ::  a3     !<
    2454     REAL(wp) ::  a4     !<
    2455     REAL(wp) ::  a5     !<
    2456     REAL(wp) ::  a6     !<
    2457     REAL(wp) ::  b0     !<
    2458     REAL(wp) ::  b1     !<
    2459     REAL(wp) ::  b10    !<
    2460     REAL(wp) ::  b11    !<
    2461     REAL(wp) ::  b2     !<
    2462     REAL(wp) ::  b20    !<
    2463     REAL(wp) ::  b21    !<
    2464     REAL(wp) ::  b3     !<
    2465     REAL(wp) ::  b4     !<
    2466     REAL(wp) ::  b5     !<
    2467     REAL(wp) ::  b6     !<
    2468     REAL(wp) ::  c1     !<
    2469     REAL(wp) ::  c2     !<
    2470     REAL(wp) ::  c3     !<
    2471     REAL(wp) ::  c4     !<
    2472     REAL(wp) ::  c5     !<
    2473     REAL(wp) ::  qrt5   !<
    2474     REAL(wp) ::  s1     !<
    2475     REAL(wp) ::  s2     !<
    2476     REAL(wp) ::  s3     !<
    2477     REAL(wp) ::  s4     !<
    2478     REAL(wp) ::  s5     !<
    2479     REAL(wp) ::  sin36  !<
    2480     REAL(wp) ::  sin45  !<
    2481     REAL(wp) ::  sin60  !<
    2482     REAL(wp) ::  sin72  !<
    2483     REAL(wp) ::  z      !<
    2484     REAL(wp) ::  zqrt5  !<
    2485     REAL(wp) ::  zsin36 !<
    2486     REAL(wp) ::  zsin45 !<
    2487     REAL(wp) ::  zsin60 !<
    2488     REAL(wp) ::  zsin72 !<
    2489 
    2490     INTEGER(iwp) ::  i     !<
    2491     INTEGER(iwp) ::  ia    !<
    2492     INTEGER(iwp) ::  ib    !<
    2493     INTEGER(iwp) ::  ibase !<
    2494     INTEGER(iwp) ::  ic    !<
    2495     INTEGER(iwp) ::  id    !<
    2496     INTEGER(iwp) ::  ie    !<
    2497     INTEGER(iwp) ::  if    !<
    2498     INTEGER(iwp) ::  ig    !<
    2499     INTEGER(iwp) ::  igo   !<
    2500     INTEGER(iwp) ::  ih    !<
    2501     INTEGER(iwp) ::  iink  !<
    2502     INTEGER(iwp) ::  ijump !<
    2503     INTEGER(iwp) ::  ipl   !<  loop index parallel loop
    2504     INTEGER(iwp) ::  j     !<
    2505     INTEGER(iwp) ::  ja    !<
    2506     INTEGER(iwp) ::  jb    !<
    2507     INTEGER(iwp) ::  jbase !<
    2508     INTEGER(iwp) ::  jc    !<
    2509     INTEGER(iwp) ::  jd    !<
    2510     INTEGER(iwp) ::  je    !<
    2511     INTEGER(iwp) ::  jf    !<
    2512     INTEGER(iwp) ::  jink  !<
    2513     INTEGER(iwp) ::  k     !<
    2514     INTEGER(iwp) ::  kb    !<
    2515     INTEGER(iwp) ::  kc    !<
    2516     INTEGER(iwp) ::  kd    !<
    2517     INTEGER(iwp) ::  ke    !<
    2518     INTEGER(iwp) ::  kf    !<
    2519     INTEGER(iwp) ::  kstop !<
    2520     INTEGER(iwp) ::  l     !<
    2521     INTEGER(iwp) ::  m     !<
     2523    REAL(wp),DIMENSION(:,:) ::  a  !<
     2524    REAL(wp),DIMENSION(:,:) ::  b  !<
     2525    REAL(wp),DIMENSION(:,:) ::  c  !<
     2526    REAL(wp),DIMENSION(:,:) ::  d  !<
     2527
     2528    REAL(wp),DIMENSION(:),INTENT(IN) ::  trigs  !<
     2529
     2530    REAL(wp) ::  a0      !<
     2531    REAL(wp) ::  a1      !<
     2532    REAL(wp) ::  a10     !<
     2533    REAL(wp) ::  a11     !<
     2534    REAL(wp) ::  a2      !<
     2535    REAL(wp) ::  a20     !<
     2536    REAL(wp) ::  a21     !<
     2537    REAL(wp) ::  a3      !<
     2538    REAL(wp) ::  a4      !<
     2539    REAL(wp) ::  a5      !<
     2540    REAL(wp) ::  a6      !<
     2541    REAL(wp) ::  b0      !<
     2542    REAL(wp) ::  b1      !<
     2543    REAL(wp) ::  b10     !<
     2544    REAL(wp) ::  b11     !<
     2545    REAL(wp) ::  b2      !<
     2546    REAL(wp) ::  b20     !<
     2547    REAL(wp) ::  b21     !<
     2548    REAL(wp) ::  b3      !<
     2549    REAL(wp) ::  b4      !<
     2550    REAL(wp) ::  b5      !<
     2551    REAL(wp) ::  b6      !<
     2552    REAL(wp) ::  c1      !<
     2553    REAL(wp) ::  c2      !<
     2554    REAL(wp) ::  c3      !<
     2555    REAL(wp) ::  c4      !<
     2556    REAL(wp) ::  c5      !<
     2557    REAL(wp) ::  qrt5    !<
     2558    REAL(wp) ::  s1      !<
     2559    REAL(wp) ::  s2      !<
     2560    REAL(wp) ::  s3      !<
     2561    REAL(wp) ::  s4      !<
     2562    REAL(wp) ::  s5      !<
     2563    REAL(wp) ::  sin36   !<
     2564    REAL(wp) ::  sin45   !<
     2565    REAL(wp) ::  sin60   !<
     2566    REAL(wp) ::  sin72   !<
     2567    REAL(wp) ::  z       !<
     2568    REAL(wp) ::  zqrt5   !<
     2569    REAL(wp) ::  zsin36  !<
     2570    REAL(wp) ::  zsin45  !<
     2571    REAL(wp) ::  zsin60  !<
     2572    REAL(wp) ::  zsin72  !<
     2573
    25222574
    25232575
     
    25312583    iink  = la * inc1
    25322584    jink  = la * inc2
    2533     ijump = (ifac-1) * iink
    2534     kstop = ( n-ifac ) / ( 2*ifac )
     2585    ijump = ( ifac - 1 ) * iink
     2586    kstop = ( n - ifac ) / ( 2 * ifac )
    25352587
    25362588    ibase = 0
     
    25532605          ib = ia + iink
    25542606          ja = 1
    2555           jb = ja + (2*m-la) * inc2
     2607          jb = ja + ( 2 * m - la ) * inc2
    25562608
    25572609          IF ( la /= m )  THEN
    25582610
    25592611             DO  l = 1, la
    2560                 i = ibase
    2561                 j = jbase
     2612                i         = ibase
     2613                j         = jbase
    25622614                c(:,ja+j) = a(:,ia+i) + a(:,ib+i)
    25632615                c(:,jb+j) = a(:,ia+i) - a(:,ib+i)
    2564                 ibase = ibase + inc1
    2565                 jbase = jbase + inc2
    2566              ENDDO
    2567 
    2568              ja = ja + jink
    2569              jink = 2 * jink
    2570              jb = jb - jink
     2616                ibase     = ibase + inc1
     2617                jbase     = jbase + inc2
     2618             ENDDO
     2619
     2620             ja    = ja + jink
     2621             jink  = 2 * jink
     2622             jb    = jb - jink
    25712623             ibase = ibase + ijump
    25722624             ijump = 2 * ijump + iink
     
    25822634                      i = ibase
    25832635                      j = jbase
    2584                       c(:,ja+j) = a(:,ia+i) + (c1*a(:,ib+i)+s1*b(:,ib+i))
    2585                       c(:,jb+j) = a(:,ia+i) - (c1*a(:,ib+i)+s1*b(:,ib+i))
    2586                       d(:,ja+j) = (c1*b(:,ib+i)-s1*a(:,ib+i)) + b(:,ia+i)
    2587                       d(:,jb+j) = (c1*b(:,ib+i)-s1*a(:,ib+i)) - b(:,ia+i)
     2636                      c(:,ja+j) = a(:,ia+i) + ( c1 * a(:,ib+i) + s1 * b(:,ib+i) )
     2637                      c(:,jb+j) = a(:,ia+i) - ( c1 * a(:,ib+i) + s1 * b(:,ib+i) )
     2638                      d(:,ja+j) = ( c1 * b(:,ib+i) - s1 * a(:,ib+i) ) + b(:,ia+i)
     2639                      d(:,jb+j) = ( c1 * b(:,ib+i) - s1 * a(:,ib+i) ) - b(:,ia+i)
    25882640                      ibase = ibase + inc1
    25892641                      jbase = jbase + inc2
     
    26112663          ELSE
    26122664
    2613              z = 1.0_wp/REAL(n,KIND=wp)
    2614              DO  l = 1, la
    2615                 i = ibase
    2616                 j = jbase
    2617                 c(:,ja+j) = z*(a(:,ia+i)+a(:,ib+i))
    2618                 c(:,jb+j) = z*(a(:,ia+i)-a(:,ib+i))
     2665             z = 1.0_wp / REAL( n, KIND = wp )
     2666             DO  l = 1, la
     2667                i = ibase
     2668                j = jbase
     2669                c(:,ja+j) = z * ( a(:,ia+i) + a(:,ib+i) )
     2670                c(:,jb+j) = z * ( a(:,ia+i) - a(:,ib+i) )
    26192671                ibase = ibase + inc1
    26202672                jbase = jbase + inc2
     
    26312683          ic = ib + iink
    26322684          ja = 1
    2633           jb = ja + (2*m-la) * inc2
     2685          jb = ja + ( 2 * m - la ) * inc2
    26342686          jc = jb
    26352687
     
    26372689
    26382690             DO  l = 1, la
    2639                 i = ibase
    2640                 j = jbase
    2641                 c(:,ja+j) = a(:,ia+i) + (a(:,ib+i)+a(:,ic+i))
    2642                 c(:,jb+j) = a(:,ia+i) - 0.5_wp*(a(:,ib+i)+a(:,ic+i))
    2643                 d(:,jb+j) = sin60*(a(:,ic+i)-a(:,ib+i))
    2644                 ibase = ibase + inc1
    2645                 jbase = jbase + inc2
    2646              ENDDO
    2647 
    2648              ja = ja + jink
    2649              jink = 2 * jink
    2650              jb = jb + jink
    2651              jc = jc - jink
     2691                i         = ibase
     2692                j         = jbase
     2693                c(:,ja+j) = a(:,ia+i) + ( a(:,ib+i) + a(:,ic+i) )
     2694                c(:,jb+j) = a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) )
     2695                d(:,jb+j) = sin60 * ( a(:,ic+i) - a(:,ib+i) )
     2696                ibase     = ibase + inc1
     2697                jbase     = jbase + inc2
     2698             ENDDO
     2699
     2700             ja    = ja + jink
     2701             jink  = 2 * jink
     2702             jb    = jb + jink
     2703             jc    = jc - jink
    26522704             ibase = ibase + ijump
    26532705             ijump = 2 * ijump + iink
     
    26672719                      i = ibase
    26682720                      j = jbase
    2669                       DO  ipl = 1, SIZE(a,1)
    2670                          a1 = (c1*a(ipl,ib+i)+s1*b(ipl,ib+i)) + (c2*a(ipl,ic+i)+s2*b(ipl,ic+i))
    2671                          b1 = (c1*b(ipl,ib+i)-s1*a(ipl,ib+i)) + (c2*b(ipl,ic+i)-s2*a(ipl,ic+i))
    2672                          a2 = a(ipl,ia+i) - 0.5_wp*a1
    2673                          b2 = b(ipl,ia+i) - 0.5_wp*b1
    2674                          a3 = sin60*((c1*a(ipl,ib+i)+s1*b(ipl,ib+i))-(c2*a(ipl,ic+i)+s2*b(ipl,ic+i)))
    2675                          b3 = sin60*((c1*b(ipl,ib+i)-s1*a(ipl,ib+i))-(c2*b(ipl,ic+i)-s2*a(ipl,ic+i)))
     2721                      DO  ipl = 1, SIZE( a, 1 )
     2722                         a1 = ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) )                              &
     2723                            + ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) )
     2724                         b1 = ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) )                              &
     2725                            + ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) )
     2726                         a2 = a(ipl,ia+i) - 0.5_wp * a1
     2727                         b2 = b(ipl,ia+i) - 0.5_wp * b1
     2728                         a3 = sin60 * ( ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) )                    &
     2729                                      - ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) ) )
     2730                         b3 = sin60 * ( ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) )                    &
     2731                                      - ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) ) )
    26762732                         c(ipl,ja+j) = a(ipl,ia+i) + a1
    26772733                         d(ipl,ja+j) = b(ipl,ia+i) + b1
     
    26792735                         d(ipl,jb+j) = b2 - a3
    26802736                         c(ipl,jc+j) = a2 - b3
    2681                          d(ipl,jc+j) = -(b2+a3)
     2737                         d(ipl,jc+j) = - ( b2 + a3 )
    26822738                      ENDDO
    26832739                      ibase = ibase + inc1
     
    26852741                   ENDDO
    26862742                   ibase = ibase + ijump
    2687                    ja = ja + jink
    2688                    jb = jb + jink
    2689                    jc = jc - jink
     2743                   ja    = ja + jink
     2744                   jb    = jb + jink
     2745                   jc    = jc - jink
    26902746                ENDDO
    26912747
     
    26982754                i = ibase
    26992755                j = jbase
    2700                 c(:,ja+j) = a(:,ia+i) + 0.5_wp*(a(:,ib+i)-a(:,ic+i))
    2701                 d(:,ja+j) = -sin60*(a(:,ib+i)+a(:,ic+i))
    2702                 c(:,jb+j) = a(:,ia+i) - (a(:,ib+i)-a(:,ic+i))
     2756                c(:,ja+j) = a(:,ia+i) + 0.5_wp * ( a(:,ib+i) - a(:,ic+i) )
     2757                d(:,ja+j) = - sin60 * ( a(:,ib+i) + a(:,ic+i) )
     2758                c(:,jb+j) = a(:,ia+i) - ( a(:,ib+i) - a(:,ic+i) )
    27032759                ibase = ibase + inc1
    27042760                jbase = jbase + inc2
     
    27072763          ELSE
    27082764
    2709              z = 1.0_wp / REAL( n, KIND=wp )
    2710              zsin60 = z*sin60
    2711              DO  l = 1, la
    2712                 i = ibase
    2713                 j = jbase
    2714                 c(:,ja+j) = z*(a(:,ia+i)+(a(:,ib+i)+a(:,ic+i)))
    2715                 c(:,jb+j) = z*(a(:,ia+i)-0.5_wp*(a(:,ib+i)+a(:,ic+i)))
    2716                 d(:,jb+j) = zsin60*(a(:,ic+i)-a(:,ib+i))
     2765             z = 1.0_wp / REAL( n, KIND = wp )
     2766             zsin60 = z * sin60
     2767             DO  l = 1, la
     2768                i = ibase
     2769                j = jbase
     2770                c(:,ja+j) = z * ( a(:,ia+i) + ( a(:,ib+i) + a(:,ic+i) ) )
     2771                c(:,jb+j) = z * ( a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) )
     2772                d(:,jb+j) = zsin60 * ( a(:,ic+i) - a(:,ib+i) )
    27172773                ibase = ibase + inc1
    27182774                jbase = jbase + inc2
     
    27302786          id = ic + iink
    27312787          ja = 1
    2732           jb = ja + (2*m-la) * inc2
    2733           jc = jb + 2*m*inc2
     2788          jb = ja + ( 2 * m - la ) * inc2
     2789          jc = jb + 2 * m * inc2
    27342790          jd = jb
    27352791
     
    27372793
    27382794             DO  l = 1, la
    2739 
    2740                 i = ibase
    2741                 j = jbase
    2742                 c(:,ja+j) = (a(:,ia+i)+a(:,ic+i)) + (a(:,ib+i)+a(:,id+i))
    2743                 c(:,jc+j) = (a(:,ia+i)+a(:,ic+i)) - (a(:,ib+i)+a(:,id+i))
     2795                i = ibase
     2796                j = jbase
     2797                c(:,ja+j) = ( a(:,ia+i) + a(:,ic+i) ) + ( a(:,ib+i) + a(:,id+i) )
     2798                c(:,jc+j) = ( a(:,ia+i) + a(:,ic+i) ) - ( a(:,ib+i) + a(:,id+i) )
    27442799                c(:,jb+j) = a(:,ia+i) - a(:,ic+i)
    27452800                d(:,jb+j) = a(:,id+i) - a(:,ib+i)
     
    27482803             ENDDO
    27492804
    2750              ja = ja + jink
    2751              jink = 2 * jink
    2752              jb = jb + jink
    2753              jc = jc - jink
    2754              jd = jd - jink
     2805             ja    = ja + jink
     2806             jink  = 2 * jink
     2807             jb    = jb + jink
     2808             jc    = jc - jink
     2809             jd    = jd - jink
    27552810             ibase = ibase + ijump
    27562811             ijump = 2 * ijump + iink
     
    27722827                      i = ibase
    27732828                      j = jbase
    2774                       DO  ipl = 1, SIZE(a,1)
    2775                          a0 = a(ipl,ia+i) + (c2*a(ipl,ic+i)+s2*b(ipl,ic+i))
    2776                          a2 = a(ipl,ia+i) - (c2*a(ipl,ic+i)+s2*b(ipl,ic+i))
    2777                          a1 = (c1*a(ipl,ib+i)+s1*b(ipl,ib+i)) + (c3*a(ipl,id+i)+s3*b(ipl,id+i))
    2778                          a3 = (c1*a(ipl,ib+i)+s1*b(ipl,ib+i)) - (c3*a(ipl,id+i)+s3*b(ipl,id+i))
    2779                          b0 = b(ipl,ia+i) + (c2*b(ipl,ic+i)-s2*a(ipl,ic+i))
    2780                          b2 = b(ipl,ia+i) - (c2*b(ipl,ic+i)-s2*a(ipl,ic+i))
    2781                          b1 = (c1*b(ipl,ib+i)-s1*a(ipl,ib+i)) + (c3*b(ipl,id+i)-s3*a(ipl,id+i))
    2782                          b3 = (c1*b(ipl,ib+i)-s1*a(ipl,ib+i)) - (c3*b(ipl,id+i)-s3*a(ipl,id+i))
     2829                      DO  ipl = 1, SIZE( a, 1 )
     2830                         a0 = a(ipl,ia+i) + ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) )
     2831                         a2 = a(ipl,ia+i) - ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) )
     2832                         a1 = ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) )                              &
     2833                            + ( c3 * a(ipl,id+i) + s3 * b(ipl,id+i) )
     2834                         a3 = ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) )                              &
     2835                            - ( c3 * a(ipl,id+i) + s3 * b(ipl,id+i) )
     2836                         b0 = b(ipl,ia+i) + ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) )
     2837                         b2 = b(ipl,ia+i) - ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) )
     2838                         b1 = ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) )                              &
     2839                            + ( c3 * b(ipl,id+i) - s3 * a(ipl,id+i) )
     2840                         b3 = ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) )                              &
     2841                            - ( c3 * b(ipl,id+i) - s3 * a(ipl,id+i) )
    27832842                         c(ipl,ja+j) = a0 + a1
    27842843                         c(ipl,jc+j) = a0 - a1
     
    27882847                         c(ipl,jd+j) = a2 - b3
    27892848                         d(ipl,jb+j) = b2 - a3
    2790                          d(ipl,jd+j) = -(b2+a3)
     2849                         d(ipl,jd+j) = - ( b2 + a3 )
    27912850                      ENDDO
    27922851                      ibase = ibase + inc1
     
    28092868                i = ibase
    28102869                j = jbase
    2811                 c(:,ja+j) = a(:,ia+i) + sin45*(a(:,ib+i)-a(:,id+i))
    2812                 c(:,jb+j) = a(:,ia+i) - sin45*(a(:,ib+i)-a(:,id+i))
    2813                 d(:,ja+j) = -a(:,ic+i) - sin45*(a(:,ib+i)+a(:,id+i))
    2814                 d(:,jb+j) = a(:,ic+i) - sin45*(a(:,ib+i)+a(:,id+i))
     2870                c(:,ja+j) =   a(:,ia+i) + sin45 * ( a(:,ib+i) - a(:,id+i) )
     2871                c(:,jb+j) =   a(:,ia+i) - sin45 * ( a(:,ib+i) - a(:,id+i) )
     2872                d(:,ja+j) = - a(:,ic+i) - sin45 * ( a(:,ib+i) + a(:,id+i) )
     2873                d(:,jb+j) =   a(:,ic+i) - sin45 * ( a(:,ib+i) + a(:,id+i) )
    28152874                ibase = ibase + inc1
    28162875                jbase = jbase + inc2
     
    28192878          ELSE
    28202879
    2821              z = 1.0_wp / REAL( n, KIND=wp )
    2822              DO  l = 1, la
    2823                 i = ibase
    2824                 j = jbase
    2825                 c(:,ja+j) = z*((a(:,ia+i)+a(:,ic+i))+(a(:,ib+i)+a(:,id+i)))
    2826                 c(:,jc+j) = z*((a(:,ia+i)+a(:,ic+i))-(a(:,ib+i)+a(:,id+i)))
    2827                 c(:,jb+j) = z*(a(:,ia+i)-a(:,ic+i))
    2828                 d(:,jb+j) = z*(a(:,id+i)-a(:,ib+i))
     2880             z = 1.0_wp / REAL( n, KIND = wp )
     2881             DO  l = 1, la
     2882                i = ibase
     2883                j = jbase
     2884                c(:,ja+j) = z * ( ( a(:,ia+i) + a(:,ic+i) ) + ( a(:,ib+i) + a(:,id+i) ) )
     2885                c(:,jc+j) = z * ( ( a(:,ia+i) + a(:,ic+i) ) - ( a(:,ib+i) + a(:,id+i) ) )
     2886                c(:,jb+j) = z * ( a(:,ia+i) - a(:,ic+i) )
     2887                d(:,jb+j) = z * ( a(:,id+i) - a(:,ib+i) )
    28292888                ibase = ibase + inc1
    28302889                jbase = jbase + inc2
     
    28432902          ie = id + iink
    28442903          ja = 1
    2845           jb = ja + (2*m-la) * inc2
    2846           jc = jb + 2*m*inc2
     2904          jb = ja + ( 2 * m - la ) * inc2
     2905          jc = jb + 2 * m * inc2
    28472906          jd = jc
    28482907          je = jb
     
    28532912                i = ibase
    28542913                j = jbase
    2855                 DO  ipl = 1, SIZE(a,1)
     2914                DO  ipl = 1, SIZE( a, 1 )
    28562915                   a1 = a(ipl,ib+i) + a(ipl,ie+i)
    28572916                   a3 = a(ipl,ib+i) - a(ipl,ie+i)
    28582917                   a2 = a(ipl,ic+i) + a(ipl,id+i)
    28592918                   a4 = a(ipl,ic+i) - a(ipl,id+i)
    2860                    a5 = a(ipl,ia+i) - 0.25_wp*(a1+a2)
    2861                    a6 = qrt5*(a1-a2)
    2862                    c(ipl,ja+j) = a(ipl,ia+i) + (a1+a2)
     2919                   a5 = a(ipl,ia+i) - 0.25_wp * ( a1 + a2 )
     2920                   a6 = qrt5 * ( a1 - a2 )
     2921                   c(ipl,ja+j) = a(ipl,ia+i) + ( a1 + a2 )
    28632922                   c(ipl,jb+j) = a5 + a6
    28642923                   c(ipl,jc+j) = a5 - a6
    2865                    d(ipl,jb+j) = -sin72*a3 - sin36*a4
    2866                    d(ipl,jc+j) = -sin36*a3 + sin72*a4
    2867                 ENDDO
    2868                 ibase = ibase + inc1
    2869                 jbase = jbase + inc2
    2870              ENDDO
    2871 
    2872              ja = ja + jink
    2873              jink = 2 * jink
    2874              jb = jb + jink
    2875              jc = jc + jink
    2876              jd = jd - jink
    2877              je = je - jink
     2924                   d(ipl,jb+j) = - sin72 * a3 - sin36 * a4
     2925                   d(ipl,jc+j) = - sin36 * a3 + sin72 * a4
     2926                ENDDO
     2927                ibase = ibase + inc1
     2928                jbase = jbase + inc2
     2929             ENDDO
     2930
     2931             ja    = ja + jink
     2932             jink  = 2 * jink
     2933             jb    = jb + jink
     2934             jc    = jc + jink
     2935             jd    = jd - jink
     2936             je    = je - jink
    28782937             ibase = ibase + ijump
    28792938             ijump = 2 * ijump + iink
     
    28992958                      j = jbase
    29002959                      DO  ipl = 1, SIZE(a,1)
    2901                          a1 = (c1*a(ipl,ib+i)+s1*b(ipl,ib+i)) + (c4*a(ipl,ie+i)+s4*b(ipl,ie+i))
    2902                          a3 = (c1*a(ipl,ib+i)+s1*b(ipl,ib+i)) - (c4*a(ipl,ie+i)+s4*b(ipl,ie+i))
    2903                          a2 = (c2*a(ipl,ic+i)+s2*b(ipl,ic+i)) + (c3*a(ipl,id+i)+s3*b(ipl,id+i))
    2904                          a4 = (c2*a(ipl,ic+i)+s2*b(ipl,ic+i)) - (c3*a(ipl,id+i)+s3*b(ipl,id+i))
    2905                          b1 = (c1*b(ipl,ib+i)-s1*a(ipl,ib+i)) + (c4*b(ipl,ie+i)-s4*a(ipl,ie+i))
    2906                          b3 = (c1*b(ipl,ib+i)-s1*a(ipl,ib+i)) - (c4*b(ipl,ie+i)-s4*a(ipl,ie+i))
    2907                          b2 = (c2*b(ipl,ic+i)-s2*a(ipl,ic+i)) + (c3*b(ipl,id+i)-s3*a(ipl,id+i))
    2908                          b4 = (c2*b(ipl,ic+i)-s2*a(ipl,ic+i)) - (c3*b(ipl,id+i)-s3*a(ipl,id+i))
    2909                          a5 = a(ipl,ia+i) - 0.25_wp*(a1+a2)
    2910                          a6 = qrt5*(a1-a2)
    2911                          b5 = b(ipl,ia+i) - 0.25_wp*(b1+b2)
    2912                          b6 = qrt5*(b1-b2)
     2960                         a1 = ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) ) + ( c4 * a(ipl,ie+i) + s4    &
     2961                              * b(ipl,ie+i) )
     2962                         a3 = ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) ) - ( c4 * a(ipl,ie+i) + s4    &
     2963                              * b(ipl,ie+i) )
     2964                         a2 = ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) ) + ( c3 * a(ipl,id+i) + s3    &
     2965                              * b(ipl,id+i) )
     2966                         a4 = ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) ) - ( c3 * a(ipl,id+i) + s3    &
     2967                              * b(ipl,id+i) )
     2968
     2969                         b1 = ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) ) + ( c4 * b(ipl,ie+i) - s4    &
     2970                              * a(ipl,ie+i) )
     2971                         b3 = ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) ) - ( c4 * b(ipl,ie+i) - s4    &
     2972                              * a(ipl,ie+i) )
     2973                         b2 = ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) ) + ( c3 * b(ipl,id+i) - s3    &
     2974                              * a(ipl,id+i) )
     2975                         b4 = ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) ) - ( c3 * b(ipl,id+i) - s3    &
     2976                              * a(ipl,id+i) )
     2977
     2978                         a5 = a(ipl,ia+i) - 0.25_wp * ( a1 + a2 )
     2979                         a6 = qrt5 * ( a1 - a2 )
     2980                         b5 = b(ipl,ia+i) - 0.25_wp * ( b1 + b2 )
     2981                         b6 = qrt5 * ( b1 - b2)
     2982
    29132983                         a10 = a5 + a6
    29142984                         a20 = a5 - a6
    29152985                         b10 = b5 + b6
    29162986                         b20 = b5 - b6
    2917                          a11 = sin72*b3 + sin36*b4
    2918                          a21 = sin36*b3 - sin72*b4
    2919                          b11 = sin72*a3 + sin36*a4
    2920                          b21 = sin36*a3 - sin72*a4
    2921                          c(ipl,ja+j) = a(ipl,ia+i) + (a1+a2)
     2987                         a11 = sin72 * b3 + sin36 * b4
     2988                         a21 = sin36 * b3 - sin72 * b4
     2989
     2990                         b11 = sin72 * a3 + sin36 * a4
     2991                         b21 = sin36 * a3 - sin72 * a4
     2992
     2993                         c(ipl,ja+j) = a(ipl,ia+i) + ( a1 + a2 )
    29222994                         c(ipl,jb+j) = a10 + a11
    29232995                         c(ipl,je+j) = a10 - a11
    29242996                         c(ipl,jc+j) = a20 + a21
    29252997                         c(ipl,jd+j) = a20 - a21
    2926                          d(ipl,ja+j) = b(ipl,ia+i) + (b1+b2)
     2998
     2999                         d(ipl,ja+j) = b(ipl,ia+i) + ( b1 + b2 )
    29273000                         d(ipl,jb+j) = b10 - b11
    2928                          d(ipl,je+j) = -(b10+b11)
     3001                         d(ipl,je+j) = - ( b10 + b11 )
    29293002                         d(ipl,jc+j) = b20 - b21
    2930                          d(ipl,jd+j) = -(b20+b21)
     3003                         d(ipl,jd+j) = - ( b20 + b21 )
    29313004                      ENDDO
    29323005                      ibase = ibase + inc1
     
    29343007                   ENDDO
    29353008                   ibase = ibase + ijump
    2936                    ja = ja + jink
    2937                    jb = jb + jink
    2938                    jc = jc + jink
    2939                    jd = jd - jink
    2940                    je = je - jink
     3009                   ja    = ja + jink
     3010                   jb    = jb + jink
     3011                   jc    = jc + jink
     3012                   jd    = jd - jink
     3013                   je    = je - jink
    29413014                ENDDO
    29423015
     
    29493022                i = ibase
    29503023                j = jbase
    2951                 DO  ipl = 1, SIZE(a,1)
     3024                DO  ipl = 1, SIZE( a, 1 )
    29523025                   a1 = a(ipl,ib+i) + a(ipl,ie+i)
    29533026                   a3 = a(ipl,ib+i) - a(ipl,ie+i)
    29543027                   a2 = a(ipl,ic+i) + a(ipl,id+i)
    29553028                   a4 = a(ipl,ic+i) - a(ipl,id+i)
    2956                    a5 = a(ipl,ia+i) + 0.25_wp*(a3-a4)
    2957                    a6 = qrt5*(a3+a4)
     3029                   a5 = a(ipl,ia+i) + 0.25_wp * ( a3 - a4 )
     3030                   a6 = qrt5 * ( a3 + a4 )
     3031
    29583032                   c(ipl,ja+j) = a5 + a6
    29593033                   c(ipl,jb+j) = a5 - a6
    2960                    c(ipl,jc+j) = a(ipl,ia+i) - (a3-a4)
    2961                    d(ipl,ja+j) = -sin36*a1 - sin72*a2
    2962                    d(ipl,jb+j) = -sin72*a1 + sin36*a2
     3034                   c(ipl,jc+j) = a(ipl,ia+i) - ( a3 - a4 )
     3035
     3036                   d(ipl,ja+j) = - sin36 * a1 - sin72 * a2
     3037                   d(ipl,jb+j) = - sin72 * a1 + sin36 * a2
    29633038                ENDDO
    29643039                ibase = ibase + inc1
     
    29683043          ELSE
    29693044
    2970              z = 1.0_wp / REAL( n, KIND=wp )
     3045             z = 1.0_wp / REAL( n, KIND = wp )
    29713046             zqrt5  = z * qrt5
    29723047             zsin36 = z * sin36
     
    29753050                i = ibase
    29763051                j = jbase
    2977                 DO  ipl = 1, SIZE(a,1)
     3052                DO  ipl = 1, SIZE( a, 1 )
    29783053                   a1 = a(ipl,ib+i) + a(ipl,ie+i)
    29793054                   a3 = a(ipl,ib+i) - a(ipl,ie+i)
    29803055                   a2 = a(ipl,ic+i) + a(ipl,id+i)
    29813056                   a4 = a(ipl,ic+i) - a(ipl,id+i)
    2982                    a5 = z*(a(ipl,ia+i)-0.25_wp*(a1+a2))
    2983                    a6 = zqrt5*(a1-a2)
    2984                    c(ipl,ja+j) = z*(a(ipl,ia+i)+(a1+a2))
     3057                   a5 = z * ( a(ipl,ia+i) - 0.25_wp * ( a1 + a2 ) )
     3058                   a6 = zqrt5 * ( a1 - a2 )
     3059
     3060                   c(ipl,ja+j) = z * ( a(ipl,ia+i) + ( a1 + a2 ) )
    29853061                   c(ipl,jb+j) = a5 + a6
    29863062                   c(ipl,jc+j) = a5 - a6
    2987                    d(ipl,jb+j) = -zsin72*a3 - zsin36*a4
    2988                    d(ipl,jc+j) = -zsin36*a3 + zsin72*a4
     3063
     3064                   d(ipl,jb+j) = - zsin72 * a3 - zsin36 * a4
     3065                   d(ipl,jc+j) = - zsin36 * a3 + zsin72 * a4
    29893066                ENDDO
    29903067                ibase = ibase + inc1
     
    30053082          if = ie + iink
    30063083          ja = 1
    3007           jb = ja + (2*m-la) * inc2
    3008           jc = jb + 2*m*inc2
    3009           jd = jc + 2*m*inc2
     3084          jb = ja + ( 2 * m - la ) * inc2
     3085          jc = jb + 2 * m * inc2
     3086          jd = jc + 2 * m * inc2
    30103087          je = jc
    30113088          jf = jb
     
    30163093                i = ibase
    30173094                j = jbase
    3018                 DO  ipl = 1, SIZE(a,1)
    3019                    a11 = (a(ipl,ic+i)+a(ipl,if+i)) + (a(ipl,ib+i)+a(ipl,ie+i))
    3020                    c(ipl,ja+j) = (a(ipl,ia+i)+a(ipl,id+i)) + a11
    3021                    c(ipl,jc+j) = (a(ipl,ia+i)+a(ipl,id+i)-0.5_wp*a11)
    3022                    d(ipl,jc+j) = sin60*((a(ipl,ic+i)+a(ipl,if+i))-(a(ipl,ib+i)+a(ipl,ie+i)))
    3023                    a11 = (a(ipl,ic+i)-a(ipl,if+i)) + (a(ipl,ie+i)-a(ipl,ib+i))
    3024                    c(ipl,jb+j) = (a(ipl,ia+i)-a(ipl,id+i)) - 0.5_wp*a11
    3025                    d(ipl,jb+j) = sin60*((a(ipl,ie+i)-a(ipl,ib+i))-(a(ipl,ic+i)-a(ipl,if+i)))
    3026                    c(ipl,jd+j) = (a(ipl,ia+i)-a(ipl,id+i)) + a11
     3095                DO  ipl = 1, SIZE( a, 1 )
     3096                   a11         = ( a(ipl,ic+i) + a(ipl,if+i) ) + ( a(ipl,ib+i) + a(ipl,ie+i) )
     3097                   c(ipl,ja+j) = ( a(ipl,ia+i) + a(ipl,id+i) ) + a11
     3098                   c(ipl,jc+j) = ( a(ipl,ia+i) + a(ipl,id+i) - 0.5_wp * a11 )
     3099                   d(ipl,jc+j) = sin60 * ( ( a(ipl,ic+i) + a(ipl,if+i) ) -                         &
     3100                                 ( a(ipl,ib+i) + a(ipl,ie+i) ) )
     3101                   a11         = ( a(ipl,ic+i) - a(ipl,if+i) ) + ( a(ipl,ie+i) - a(ipl,ib+i) )
     3102                   c(ipl,jb+j) = ( a(ipl,ia+i) - a(ipl,id+i) ) - 0.5_wp * a11
     3103                   d(ipl,jb+j) = sin60 * ( ( a(ipl,ie+i) - a(ipl,ib+i) ) -                         &
     3104                                 ( a(ipl,ic+i) - a(ipl,if+i) ) )
     3105                   c(ipl,jd+j) = ( a(ipl,ia+i) - a(ipl,id+i) ) + a11
    30273106                END DO
    30283107                ibase = ibase + inc1
     
    30613140                      i = ibase
    30623141                      j = jbase
    3063                       DO  ipl = 1, SIZE(a,1)
    3064                          a1 = c1*a(ipl,ib+i) + s1*b(ipl,ib+i)
    3065                          b1 = c1*b(ipl,ib+i) - s1*a(ipl,ib+i)
    3066                          a2 = c2*a(ipl,ic+i) + s2*b(ipl,ic+i)
    3067                          b2 = c2*b(ipl,ic+i) - s2*a(ipl,ic+i)
    3068                          a3 = c3*a(ipl,id+i) + s3*b(ipl,id+i)
    3069                          b3 = c3*b(ipl,id+i) - s3*a(ipl,id+i)
    3070                          a4 = c4*a(ipl,ie+i) + s4*b(ipl,ie+i)
    3071                          b4 = c4*b(ipl,ie+i) - s4*a(ipl,ie+i)
    3072                          a5 = c5*a(ipl,if+i) + s5*b(ipl,if+i)
    3073                          b5 = c5*b(ipl,if+i) - s5*a(ipl,if+i)
    3074                          a11 = (a2+a5) + (a1+a4)
    3075                          a20 = (a(ipl,ia+i)+a3) - 0.5_wp*a11
    3076                          a21 = sin60*((a2+a5)-(a1+a4))
    3077                          b11 = (b2+b5) + (b1+b4)
    3078                          b20 = (b(ipl,ia+i)+b3) - 0.5_wp*b11
    3079                          b21 = sin60*((b2+b5)-(b1+b4))
    3080                          c(ipl,ja+j) = (a(ipl,ia+i)+a3) + a11
    3081                          d(ipl,ja+j) = (b(ipl,ia+i)+b3) + b11
     3142                      DO  ipl = 1, SIZE( a, 1 )
     3143                         a1 = c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i)
     3144                         b1 = c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i)
     3145                         a2 = c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i)
     3146                         b2 = c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i)
     3147                         a3 = c3 * a(ipl,id+i) + s3 * b(ipl,id+i)
     3148                         b3 = c3 * b(ipl,id+i) - s3 * a(ipl,id+i)
     3149                         a4 = c4 * a(ipl,ie+i) + s4 * b(ipl,ie+i)
     3150                         b4 = c4 * b(ipl,ie+i) - s4 * a(ipl,ie+i)
     3151                         a5 = c5 * a(ipl,if+i) + s5 * b(ipl,if+i)
     3152                         b5 = c5 * b(ipl,if+i) - s5 * a(ipl,if+i)
     3153                         a11 = ( a2 + a5 ) + ( a1 + a4 )
     3154                         a20 = ( a(ipl,ia+i) + a3 ) - 0.5_wp * a11
     3155                         a21 = sin60 * ( ( a2 + a5 ) - ( a1 + a4 ) )
     3156                         b11 = ( b2 + b5 ) + ( b1 + b4 )
     3157                         b20 = ( b(ipl,ia+i) + b3 ) - 0.5_wp * b11
     3158                         b21 = sin60 * ( ( b2 + b5 ) - ( b1 + b4 ) )
     3159                         c(ipl,ja+j) = ( a(ipl,ia+i) + a3 ) + a11
     3160                         d(ipl,ja+j) = ( b(ipl,ia+i) + b3 ) + b11
    30823161                         c(ipl,jc+j) = a20 - b21
    30833162                         d(ipl,jc+j) = a21 + b20
    30843163                         c(ipl,je+j) = a20 + b21
    30853164                         d(ipl,je+j) = a21 - b20
    3086                          a11 = (a2-a5) + (a4-a1)
    3087                          a20 = (a(ipl,ia+i)-a3) - 0.5_wp*a11
    3088                          a21 = sin60*((a4-a1)-(a2-a5))
    3089                          b11 = (b5-b2) - (b4-b1)
    3090                          b20 = (b3-b(ipl,ia+i)) - 0.5_wp*b11
    3091                          b21 = sin60*((b5-b2)+(b4-b1))
     3165                         a11 = ( a2 - a5 ) + ( a4 - a1 )
     3166                         a20 = ( a(ipl,ia+i) - a3 ) - 0.5_wp * a11
     3167                         a21 = sin60 * ( ( a4 - a1 ) - ( a2 - a5 ) )
     3168                         b11 = ( b5 - b2 ) - ( b4 - b1 )
     3169                         b20 = ( b3 - b(ipl,ia+i) ) - 0.5_wp * b11
     3170                         b21 = sin60 * ( ( b5 - b2 ) + ( b4 - b1 ) )
    30923171                         c(ipl,jb+j) = a20 - b21
    30933172                         d(ipl,jb+j) = a21 - b20
    3094                          c(ipl,jd+j) = a11 + (a(ipl,ia+i)-a3)
    3095                          d(ipl,jd+j) = b11 + (b3-b(ipl,ia+i))
     3173                         c(ipl,jd+j) = a11 + ( a(ipl,ia+i) - a3 )
     3174                         d(ipl,jd+j) = b11 + ( b3 - b(ipl,ia+i) )
    30963175                         c(ipl,jf+j) = a20 + b21
    30973176                         d(ipl,jf+j) = a21 + b20
     
    31173196                i = ibase
    31183197                j = jbase
    3119                 c(:,ja+j) = (a(:,ia+i)+0.5_wp*(a(:,ic+i)-a(:,ie+i))) + sin60*(a(:,ib+i)-a(:,if+i))
    3120                 d(:,ja+j) = -(a(:,id+i)+0.5_wp*(a(:,ib+i)+a(:,if+i))) - sin60*(a(:,ic+i)+a(:,ie+i))
    3121                 c(:,jb+j) = a(:,ia+i) - (a(:,ic+i)-a(:,ie+i))
    3122                 d(:,jb+j) = a(:,id+i) - (a(:,ib+i)+a(:,if+i))
    3123                 c(:,jc+j) = (a(:,ia+i)+0.5_wp*(a(:,ic+i)-a(:,ie+i))) - sin60*(a(:,ib+i)-a(:,if+i))
    3124                 d(:,jc+j) = -(a(:,id+i)+0.5_wp*(a(:,ib+i)+a(:,if+i))) + sin60*(a(:,ic+i)+a(:,ie+i))
     3198                c(:,ja+j) = ( a(:,ia+i) + 0.5_wp * ( a(:,ic+i) - a(:,ie+i) ) ) + sin60             &
     3199                          * ( a(:,ib+i) - a(:,if+i) )
     3200                d(:,ja+j) = - ( a(:,id+i) + 0.5_wp * ( a(:,ib+i) + a(:,if+i) ) ) - sin60           &
     3201                            * ( a(:,ic+i) + a(:,ie+i) )
     3202                c(:,jb+j) = a(:,ia+i) - ( a(:,ic+i) - a(:,ie+i) )
     3203                d(:,jb+j) = a(:,id+i) - ( a(:,ib+i) + a(:,if+i) )
     3204                c(:,jc+j) = ( a(:,ia+i) + 0.5_wp * ( a(:,ic+i) - a(:,ie+i) ) ) - sin60             &
     3205                          * ( a(:,ib+i) - a(:,if+i) )
     3206                d(:,jc+j) = - ( a(:,id+i) + 0.5_wp * ( a(:,ib+i) + a(:,if+i) ) ) + sin60           &
     3207                            * ( a(:,ic+i) + a(:,ie+i) )
    31253208                ibase = ibase + inc1
    31263209                jbase = jbase + inc2
     
    31293212          ELSE
    31303213
    3131              z = 1.0_wp/REAL(n,KIND=wp)
    3132              zsin60 = z*sin60
    3133              DO  l = 1, la
    3134                 i = ibase
    3135                 j = jbase
    3136                 DO  ipl = 1, SIZE(a,1)
    3137                    a11 = (a(ipl,ic+i)+a(ipl,if+i)) + (a(ipl,ib+i)+a(ipl,ie+i))
    3138                    c(ipl,ja+j) = z*((a(ipl,ia+i)+a(ipl,id+i))+a11)
    3139                    c(ipl,jc+j) = z*((a(ipl,ia+i)+a(ipl,id+i))-0.5_wp*a11)
    3140                    d(ipl,jc+j) = zsin60*((a(ipl,ic+i)+a(ipl,if+i))-(a(ipl,ib+i)+a(ipl,ie+i)))
    3141                    a11 = (a(ipl,ic+i)-a(ipl,if+i)) + (a(ipl,ie+i)-a(ipl,ib+i))
    3142                    c(ipl,jb+j) = z*((a(ipl,ia+i)-a(ipl,id+i))-0.5_wp*a11)
    3143                    d(ipl,jb+j) = zsin60*((a(ipl,ie+i)-a(ipl,ib+i))-(a(ipl,ic+i)-a(ipl,if+i)))
    3144                    c(ipl,jd+j) = z*((a(ipl,ia+i)-a(ipl,id+i))+a11)
     3214             z = 1.0_wp / REAL( n, KIND = wp )
     3215             zsin60 = z * sin60
     3216             DO  l = 1, la
     3217                i = ibase
     3218                j = jbase
     3219                DO  ipl = 1, SIZE( a, 1 )
     3220                   a11 = ( a (ipl,ic+i) + a(ipl,if+i) ) + ( a(ipl,ib+i) + a(ipl,ie+i) )
     3221                   c(ipl,ja+j) = z * ( ( a(ipl,ia+i) + a(ipl,id+i) ) + a11 )
     3222                   c(ipl,jc+j) = z * ( ( a(ipl,ia+i) + a(ipl,id+i) ) - 0.5_wp * a11 )
     3223                   d(ipl,jc+j) = zsin60 * ( ( a(ipl,ic+i) + a(ipl,if+i) )                          &
     3224                                          - ( a(ipl,ib+i) + a(ipl,ie+i) ) )
     3225                   a11 = ( a(ipl,ic+i) - a(ipl,if+i) )  + ( a(ipl,ie+i) - a(ipl,ib+i) )
     3226                   c(ipl,jb+j) = z * ( ( a(ipl,ia+i) - a(ipl,id+i) ) - 0.5_wp * a11 )
     3227                   d(ipl,jb+j) = zsin60 * ( ( a(ipl,ie+i) - a(ipl,ib+i) )                          &
     3228                                          - ( a(ipl,ic+i) - a(ipl,if+i) ) )
     3229                   c(ipl,jd+j) = z * ( ( a(ipl,ia+i) - a(ipl,id+i) ) + a11 )
    31453230                ENDDO
    31463231                ibase = ibase + inc1
     
    31693254          ja = 1
    31703255          jb = ja + la * inc2
    3171           jc = jb + 2*m*inc2
    3172           jd = jc + 2*m*inc2
    3173           je = jd + 2*m*inc2
    3174           z = 1.0_wp / REAL( n, KIND=wp )
     3256          jc = jb + 2 * m * inc2
     3257          jd = jc + 2 * m * inc2
     3258          je = jd + 2 * m * inc2
     3259          z = 1.0_wp / REAL( n, KIND = wp )
    31753260          zsin45 = z * SQRT( 0.5_wp )
    31763261
     
    31783263             i = ibase
    31793264             j = jbase
    3180              c(:,ja+j) = z*(((a(:,ia+i)+a(:,ie+i))+(a(:,ic+i)+a(:,ig+i)))+((a(:,id+i)+ a(:,ih+i))+(a(:,ib+i)+a(:,if+i))))
    3181              c(:,je+j) = z*(((a(:,ia+i)+a(:,ie+i))+(a(:,ic+i)+a(:,ig+i)))-((a(:,id+i)+ a(:,ih+i))+(a(:,ib+i)+a(:,if+i))))
    3182              c(:,jc+j) = z*((a(:,ia+i)+a(:,ie+i))-(a(:,ic+i)+a(:,ig+i)))
    3183              d(:,jc+j) = z*((a(:,id+i)+a(:,ih+i))-(a(:,ib+i)+a(:,if+i)))
    3184              c(:,jb+j) = z*(a(:,ia+i)-a(:,ie+i)) + zsin45*((a(:,ih+i)-a(:,id+i))-(a(:,if+i)-a(:,ib+i)))
    3185              c(:,jd+j) = z*(a(:,ia+i)-a(:,ie+i)) - zsin45*((a(:,ih+i)-a(:,id+i))-(a(:,if+i)-a(:,ib+i)))
    3186              d(:,jb+j) = zsin45*((a(:,ih+i)-a(:,id+i))+(a(:,if+i)-a(:,ib+i))) + z*(a(:,ig+i)-a(:,ic+i))
    3187              d(:,jd+j) = zsin45*((a(:,ih+i)-a(:,id+i))+(a(:,if+i)-a(:,ib+i))) - z*(a(:,ig+i)-a(:,ic+i))
     3265             c(:,ja+j) = z * ( ( ( a(:,ia+i) + a(:,ie+i) ) + ( a(:,ic+i) + a(:,ig+i) ) )           &
     3266                         + ( ( a(:,id+i) + a(:,ih+i) ) + ( a(:,ib+i) + a(:,if+i) ) ) )
     3267             c(:,je+j) = z * ( ( ( a(:,ia+i) + a(:,ie+i) ) + ( a(:,ic+i) + a(:,ig+i) ) )           &
     3268                         - ( ( a(:,id+i) + a(:,ih+i) ) + ( a(:,ib+i) + a(:,if+i) ) ) )
     3269             c(:,jc+j) = z * ( ( a(:,ia+i) + a(:,ie+i) ) - ( a(:,ic+i) + a(:,ig+i) ) )
     3270             d(:,jc+j) = z * ( ( a(:,id+i) + a(:,ih+i) ) - ( a(:,ib+i) + a(:,if+i) ) )
     3271             c(:,jb+j) = z * ( a(:,ia+i) - a(:,ie+i) ) + zsin45 * ( ( a(:,ih+i) - a(:,id+i) )      &
     3272                         - ( a(:,if+i) - a(:,ib+i) ) )
     3273             c(:,jd+j) = z * ( a(:,ia+i) - a(:,ie+i) ) - zsin45 * ( ( a(:,ih+i) - a(:,id+i) )      &
     3274                         - ( a(:,if+i) - a(:,ib+i) ) )
     3275             d(:,jb+j) = zsin45 * ( ( a(:,ih+i) - a(:,id+i) ) + ( a(:,if+i) - a(:,ib+i) ) )        &
     3276                         + z * ( a(:,ig+i) - a(:,ic+i) )
     3277             d(:,jd+j) = zsin45 * ( ( a(:,ih+i) - a(:,id+i) ) + ( a(:,if+i) - a(:,ib+i) ) )        &
     3278                         - z * ( a(:,ig+i) - a(:,ic+i) )
    31883279             ibase = ibase + inc1
    31893280             jbase = jbase + inc2
     
    31943285 END SUBROUTINE qpassm_vec
    31953286
    3196 !------------------------------------------------------------------------------!
     3287!--------------------------------------------------------------------------------------------------!
    31973288! Description:
    31983289! ------------
    31993290!> Same as qpassm, but for backward fft
    3200 !------------------------------------------------------------------------------!
     3291!--------------------------------------------------------------------------------------------------!
    32013292 SUBROUTINE rpassm_vec(a, b, c, d, trigs, inc1, inc2, n, ifac, la, ierr )
    32023293
     
    32053296    IMPLICIT NONE
    32063297
    3207     INTEGER(iwp) ::  ierr !<
    3208     INTEGER(iwp) ::  ifac !<
    3209     INTEGER(iwp) ::  inc1 !<
    3210     INTEGER(iwp) ::  inc2 !<
    3211     INTEGER(iwp) ::  la   !<
    3212     INTEGER(iwp) ::  n    !<
     3298
     3299    INTEGER(iwp) ::  i      !<
     3300    INTEGER(iwp) ::  ia     !<
     3301    INTEGER(iwp) ::  ib     !<
     3302    INTEGER(iwp) ::  ibase  !<
     3303    INTEGER(iwp) ::  ic     !<
     3304    INTEGER(iwp) ::  id     !<
     3305    INTEGER(iwp) ::  ie     !<
     3306    INTEGER(iwp) ::  ierr   !<
     3307    INTEGER(iwp) ::  ifac   !<
     3308    INTEGER(iwp) ::  if     !<
     3309    INTEGER(iwp) ::  igo    !<
     3310    INTEGER(iwp) ::  iink   !<
     3311    INTEGER(iwp) ::  inc1   !<
     3312    INTEGER(iwp) ::  inc2   !<
     3313    INTEGER(iwp) ::  ipl    !<  loop index parallel loop
     3314    INTEGER(iwp) ::  j      !<
     3315    INTEGER(iwp) ::  ja     !<
     3316    INTEGER(iwp) ::  jb     !<
     3317    INTEGER(iwp) ::  jbase  !<
     3318    INTEGER(iwp) ::  jc     !<
     3319    INTEGER(iwp) ::  jd     !<
     3320    INTEGER(iwp) ::  je     !<
     3321    INTEGER(iwp) ::  jf     !<
     3322    INTEGER(iwp) ::  jg     !<
     3323    INTEGER(iwp) ::  jh     !<
     3324    INTEGER(iwp) ::  jink   !<
     3325    INTEGER(iwp) ::  jump   !<
     3326    INTEGER(iwp) ::  k      !<
     3327    INTEGER(iwp) ::  kb     !<
     3328    INTEGER(iwp) ::  kc     !<
     3329    INTEGER(iwp) ::  kd     !<
     3330    INTEGER(iwp) ::  ke     !<
     3331    INTEGER(iwp) ::  kf     !<
     3332    INTEGER(iwp) ::  kstop  !<
     3333    INTEGER(iwp) ::  l      !<
     3334    INTEGER(iwp) ::  la     !<
     3335    INTEGER(iwp) ::  m      !<
     3336    INTEGER(iwp) ::  n      !<
    32133337!
    32143338!-- Arrays are dimensioned with n
    3215     REAL(wp),DIMENSION(:,:) ::  a     !<
    3216     REAL(wp),DIMENSION(:,:) ::  b     !<
    3217     REAL(wp),DIMENSION(:,:) ::  c     !<
    3218     REAL(wp),DIMENSION(:,:) ::  d     !<
    3219     REAL(wp),DIMENSION(:),INTENT(IN) ::  trigs !<
    3220 
    3221     REAL(wp) ::  c1     !<
    3222     REAL(wp) ::  c2     !<
    3223     REAL(wp) ::  c3     !<
    3224     REAL(wp) ::  c4     !<
    3225     REAL(wp) ::  c5     !<
    3226     REAL(wp) ::  qqrt5  !<
    3227     REAL(wp) ::  qrt5   !<
    3228     REAL(wp) ::  s1     !<
    3229     REAL(wp) ::  s2     !<
    3230     REAL(wp) ::  s3     !<
    3231     REAL(wp) ::  s4     !<
    3232     REAL(wp) ::  s5     !<
    3233     REAL(wp) ::  sin36  !<
    3234     REAL(wp) ::  sin45  !<
    3235     REAL(wp) ::  sin60  !<
    3236     REAL(wp) ::  sin72  !<
    3237     REAL(wp) ::  ssin36 !<
    3238     REAL(wp) ::  ssin45 !<
    3239     REAL(wp) ::  ssin60 !<
    3240     REAL(wp) ::  ssin72 !<
    3241 
    3242     INTEGER(iwp) ::  i     !<
    3243     INTEGER(iwp) ::  ia    !<
    3244     INTEGER(iwp) ::  ib    !<
    3245     INTEGER(iwp) ::  ibase !<
    3246     INTEGER(iwp) ::  ic    !<
    3247     INTEGER(iwp) ::  id    !<
    3248     INTEGER(iwp) ::  ie    !<
    3249     INTEGER(iwp) ::  if    !<
    3250     INTEGER(iwp) ::  igo   !<
    3251     INTEGER(iwp) ::  iink  !<
    3252     INTEGER(iwp) ::  ipl   !<  loop index parallel loop
    3253     INTEGER(iwp) ::  j     !<
    3254     INTEGER(iwp) ::  ja    !<
    3255     INTEGER(iwp) ::  jb    !<
    3256     INTEGER(iwp) ::  jbase !<
    3257     INTEGER(iwp) ::  jc    !<
    3258     INTEGER(iwp) ::  jd    !<
    3259     INTEGER(iwp) ::  je    !<
    3260     INTEGER(iwp) ::  jf    !<
    3261     INTEGER(iwp) ::  jg    !<
    3262     INTEGER(iwp) ::  jh    !<
    3263     INTEGER(iwp) ::  jink  !<
    3264     INTEGER(iwp) ::  jump  !<
    3265     INTEGER(iwp) ::  k     !<
    3266     INTEGER(iwp) ::  kb    !<
    3267     INTEGER(iwp) ::  kc    !<
    3268     INTEGER(iwp) ::  kd    !<
    3269     INTEGER(iwp) ::  ke    !<
    3270     INTEGER(iwp) ::  kf    !<
    3271     INTEGER(iwp) ::  kstop !<
    3272     INTEGER(iwp) ::  l     !<
    3273     INTEGER(iwp) ::  m     !<
    3274 
    3275     REAL(wp) ::  a10       !<
    3276     REAL(wp) ::  a11       !<
    3277     REAL(wp) ::  a20       !<
    3278     REAL(wp) ::  a21       !<
    3279     REAL(wp) ::  b10       !<
    3280     REAL(wp) ::  b11       !<
    3281     REAL(wp) ::  b20       !<
    3282     REAL(wp) ::  b21       !<
     3339    REAL(wp),DIMENSION(:,:) ::  a  !<
     3340    REAL(wp),DIMENSION(:,:) ::  b  !<
     3341    REAL(wp),DIMENSION(:,:) ::  c  !<
     3342    REAL(wp),DIMENSION(:,:) ::  d  !<
     3343    REAL(wp),DIMENSION(:),INTENT(IN) ::  trigs  !<
     3344
     3345    REAL(wp) ::  a10     !<
     3346    REAL(wp) ::  a11     !<
     3347    REAL(wp) ::  a20     !<
     3348    REAL(wp) ::  a21     !<
     3349    REAL(wp) ::  b10     !<
     3350    REAL(wp) ::  b11     !<
     3351    REAL(wp) ::  b20     !<
     3352    REAL(wp) ::  b21     !<
     3353    REAL(wp) ::  c1      !<
     3354    REAL(wp) ::  c2      !<
     3355    REAL(wp) ::  c3      !<
     3356    REAL(wp) ::  c4      !<
     3357    REAL(wp) ::  c5      !<
     3358    REAL(wp) ::  qqrt5   !<
     3359    REAL(wp) ::  qrt5    !<
     3360    REAL(wp) ::  s1      !<
     3361    REAL(wp) ::  s2      !<
     3362    REAL(wp) ::  s3      !<
     3363    REAL(wp) ::  s4      !<
     3364    REAL(wp) ::  s5      !<
     3365    REAL(wp) ::  sin36   !<
     3366    REAL(wp) ::  sin45   !<
     3367    REAL(wp) ::  sin60   !<
     3368    REAL(wp) ::  sin72   !<
     3369    REAL(wp) ::  ssin36  !<
     3370    REAL(wp) ::  ssin45  !<
     3371    REAL(wp) ::  ssin60  !<
     3372    REAL(wp) ::  ssin72  !<
     3373
    32833374
    32843375
     
    32923383    iink = la * inc1
    32933384    jink = la * inc2
    3294     jump = (ifac-1) * jink
    3295     kstop = (n-ifac) / (2*ifac)
     3385    jump = ( ifac - 1 ) * jink
     3386    kstop = ( n - ifac ) / ( 2 * ifac )
    32963387
    32973388    ibase = 0
     
    33113402
    33123403          ia = 1
    3313           ib = ia + (2*m-la) * inc1
     3404          ib = ia + ( 2 * m - la ) * inc1
    33143405          ja = 1
    33153406          jb = ja + jink
     
    33453436                      c(:,ja+j) = a(:,ia+i) + a(:,ib+i)
    33463437                      d(:,ja+j) = b(:,ia+i) - b(:,ib+i)
    3347                       c(:,jb+j) = c1*(a(:,ia+i)-a(:,ib+i)) - s1*(b(:,ia+i)+b(:,ib+i))
    3348                       d(:,jb+j) = s1*(a(:,ia+i)-a(:,ib+i)) + c1*(b(:,ia+i)+b(:,ib+i))
     3438                      c(:,jb+j) = c1 * ( a(:,ia+i) - a(:,ib+i) ) - s1 * ( b(:,ia+i) + b(:,ib+i) )
     3439                      d(:,jb+j) = s1 * ( a(:,ia+i) - a(:,ib+i) ) + c1 * ( b(:,ia+i) + b(:,ib+i) )
    33493440                      ibase = ibase + inc1
    33503441                      jbase = jbase + inc2
     
    33643455                j = jbase
    33653456                c(:,ja+j) = a(:,ia+i)
    3366                 c(:,jb+j) = -b(:,ia+i)
     3457                c(:,jb+j) = - b(:,ia+i)
    33673458                ibase = ibase + inc1
    33683459                jbase = jbase + inc2
     
    33743465                i = ibase
    33753466                j = jbase
    3376                 c(:,ja+j) = 2.0_wp*(a(:,ia+i)+a(:,ib+i))
    3377                 c(:,jb+j) = 2.0_wp*(a(:,ia+i)-a(:,ib+i))
     3467                c(:,ja+j) = 2.0_wp * ( a(:,ia+i) + a(:,ib+i) )
     3468                c(:,jb+j) = 2.0_wp * ( a(:,ia+i) - a(:,ib+i) )
    33783469                ibase = ibase + inc1
    33793470                jbase = jbase + inc2
     
    33873478
    33883479          ia = 1
    3389           ib = ia + (2*m-la) * inc1
     3480          ib = ia + ( 2 * m - la ) * inc1
    33903481          ic = ib
    33913482          ja = 1
     
    33993490                j = jbase
    34003491                c(:,ja+j) = a(:,ia+i) + a(:,ib+i)
    3401                 c(:,jb+j) = (a(:,ia+i)-0.5_wp*a(:,ib+i)) - (sin60*(b(:,ib+i)))
    3402                 c(:,jc+j) = (a(:,ia+i)-0.5_wp*a(:,ib+i)) + (sin60*(b(:,ib+i)))
     3492                c(:,jb+j) = ( a(:,ia+i) - 0.5_wp * a(:,ib+i) ) - ( sin60 * ( b(:,ib+i) ) )
     3493                c(:,jc+j) = ( a(:,ia+i) - 0.5_wp * a(:,ib+i) ) + ( sin60 * ( b(:,ib+i) ) )
    34033494                ibase = ibase + inc1
    34043495                jbase = jbase + inc2
     
    34243515                      i = ibase
    34253516                      j = jbase
    3426                       c(:,ja+j) = a(:,ia+i) + (a(:,ib+i)+a(:,ic+i))
    3427                       d(:,ja+j) = b(:,ia+i) + (b(:,ib+i)-b(:,ic+i))
    3428                       c(:,jb+j) = c1*((a(:,ia+i)-0.5_wp*(a(:,ib+i)+a(:,ic+i)))-(sin60*(b(:,ib+i)+ b(:,ic+i)))) &
    3429                                 - s1*((b(:,ia+i)-0.5_wp*(b(:,ib+i)-b(:,ic+i)))+(sin60*(a(:,ib+i)- a(:,ic+i))))
    3430                       d(:,jb+j) = s1*((a(:,ia+i)-0.5_wp*(a(:,ib+i)+a(:,ic+i)))-(sin60*(b(:,ib+i)+ b(:,ic+i)))) &
    3431                                 + c1*((b(:,ia+i)-0.5_wp*(b(:,ib+i)-b(:,ic+i)))+(sin60*(a(:,ib+i)- a(:,ic+i))))
    3432                       c(:,jc+j) = c2*((a(:,ia+i)-0.5_wp*(a(:,ib+i)+a(:,ic+i)))+(sin60*(b(:,ib+i)+ b(:,ic+i)))) &
    3433                                 - s2*((b(:,ia+i)-0.5_wp*(b(:,ib+i)-b(:,ic+i)))-(sin60*(a(:,ib+i)- a(:,ic+i))))
    3434                       d(:,jc+j) = s2*((a(:,ia+i)-0.5_wp*(a(:,ib+i)+a(:,ic+i)))+(sin60*(b(:,ib+i)+ b(:,ic+i)))) &
    3435                                 + c2*((b(:,ia+i)-0.5_wp*(b(:,ib+i)-b(:,ic+i)))-(sin60*(a(:,ib+i)- a(:,ic+i))))
     3517                      c(:,ja+j) = a(:,ia+i) + ( a(:,ib+i) + a(:,ic+i) )
     3518                      d(:,ja+j) = b(:,ia+i) + ( b(:,ib+i) - b(:,ic+i) )
     3519                      c(:,jb+j) = c1 * ( ( a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) )        &
     3520                                  - ( sin60 * ( b(:,ib+i) + b(:,ic+i) ) ) ) - s1 * ( ( b(:,ia+i)   &
     3521                                  - 0.5_wp * ( b(:,ib+i) - b(:,ic+i) ) ) + ( sin60 * ( a(:,ib+i)   &
     3522                                  - a(:,ic+i) ) ) )
     3523                      d(:,jb+j) = s1 * ( ( a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) )        &
     3524                                  - ( sin60 * ( b(:,ib+i) + b(:,ic+i) ) ) ) + c1 * ( ( b(:,ia+i)   &
     3525                                  - 0.5_wp * ( b(:,ib+i) - b(:,ic+i) ) ) + ( sin60 * ( a(:,ib+i)   &
     3526                                  - a(:,ic+i) ) ) )
     3527                      c(:,jc+j) = c2 * ( ( a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) )        &
     3528                                  + ( sin60 * ( b(:,ib+i) + b(:,ic+i) ) ) ) - s2 * ( ( b(:,ia+i)   &
     3529                                  - 0.5_wp * ( b(:,ib+i) - b(:,ic+i) ) ) - ( sin60 * ( a(:,ib+i)   &
     3530                                  - a(:,ic+i) ) ) )
     3531                      d(:,jc+j) = s2 * ( ( a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) )        &
     3532                                  + ( sin60 * ( b(:,ib+i) + b(:,ic+i) ) ) ) + c2 * ( ( b(:,ia+i)   &
     3533                                  - 0.5_wp * ( b(:,ib+i) - b(:,ic+i) ) ) - ( sin60 * ( a(:,ib+i)   &
     3534                                  - a(:,ic+i) ) ) )
    34363535                      ibase = ibase + inc1
    34373536                      jbase = jbase + inc2
     
    34523551                j = jbase
    34533552                c(:,ja+j) = a(:,ia+i) + a(:,ib+i)
    3454                 c(:,jb+j) = (0.5_wp*a(:,ia+i)-a(:,ib+i)) - (sin60*b(:,ia+i))
    3455                 c(:,jc+j) = -(0.5_wp*a(:,ia+i)-a(:,ib+i)) - (sin60*b(:,ia+i))
     3553                c(:,jb+j) = ( 0.5_wp * a(:,ia+i) - a(:,ib+i) ) - ( sin60 * b(:,ia+i) )
     3554                c(:,jc+j) = - ( 0.5_wp * a(:,ia+i) - a(:,ib+i) ) - ( sin60 * b(:,ia+i) )
    34563555                ibase = ibase + inc1
    34573556                jbase = jbase + inc2
     
    34643563                i = ibase
    34653564                j = jbase
    3466                 c(:,ja+j) = 2.0_wp*(a(:,ia+i)+a(:,ib+i))
    3467                 c(:,jb+j) = (2.0_wp*a(:,ia+i)-a(:,ib+i)) - (ssin60*b(:,ib+i))
    3468                 c(:,jc+j) = (2.0_wp*a(:,ia+i)-a(:,ib+i)) + (ssin60*b(:,ib+i))
     3565                c(:,ja+j) = 2.0_wp * ( a(:,ia+i) + a(:,ib+i) )
     3566                c(:,jb+j) = ( 2.0_wp * a(:,ia+i) - a(:,ib+i) ) - ( ssin60 * b(:,ib+i) )
     3567                c(:,jc+j) = ( 2.0_wp * a(:,ia+i) - a(:,ib+i) ) + ( ssin60 * b(:,ib+i) )
    34693568                ibase = ibase + inc1
    34703569                jbase = jbase + inc2
     
    34783577
    34793578          ia = 1
    3480           ib = ia + (2*m-la) * inc1
    3481           ic = ib + 2*m*inc1
     3579          ib = ia + ( 2 * m - la ) * inc1
     3580          ic = ib + 2 * m * inc1
    34823581          id = ib
    34833582          ja = 1
     
    34913590                i = ibase
    34923591                j = jbase
    3493                 c(:,ja+j) = (a(:,ia+i)+a(:,ic+i)) + a(:,ib+i)
    3494                 c(:,jb+j) = (a(:,ia+i)-a(:,ic+i)) - b(:,ib+i)
    3495                 c(:,jc+j) = (a(:,ia+i)+a(:,ic+i)) - a(:,ib+i)
    3496                 c(:,jd+j) = (a(:,ia+i)-a(:,ic+i)) + b(:,ib+i)
     3592                c(:,ja+j) = ( a(:,ia+i) + a(:,ic+i) ) + a(:,ib+i)
     3593                c(:,jb+j) = ( a(:,ia+i) - a(:,ic+i) ) - b(:,ib+i)
     3594                c(:,jc+j) = ( a(:,ia+i) + a(:,ic+i) ) - a(:,ib+i)
     3595                c(:,jd+j) = ( a(:,ia+i) - a(:,ic+i) ) + b(:,ib+i)
    34973596                ibase = ibase + inc1
    34983597                jbase = jbase + inc2
     
    35223621                      i = ibase
    35233622                      j = jbase
    3524                       c(:,ja+j) = (a(:,ia+i)+a(:,ic+i)) + (a(:,ib+i)+a(:,id+i))
    3525                       d(:,ja+j) = (b(:,ia+i)-b(:,ic+i)) + (b(:,ib+i)-b(:,id+i))
    3526                       c(:,jc+j) = c2*((a(:,ia+i)+a(:,ic+i))-(a(:,ib+i)+a(:,id+i))) - s2*((b(:,ia+i)-b(:,ic+i))&
    3527                                 -(b(:,ib+i)-b(:,id+i)))
    3528                       d(:,jc+j) = s2*((a(:,ia+i)+a(:,ic+i))-(a(:,ib+i)+a(:,id+i))) + c2*((b(:,ia+i)-b(:,ic+i))&
    3529                                 -(b(:,ib+i)-b(:,id+i)))
    3530                       c(:,jb+j) = c1*((a(:,ia+i)-a(:,ic+i))-(b(:,ib+i)+b(:,id+i))) - s1*((b(:,ia+i)+b(:,ic+i))&
    3531                                 +(a(:,ib+i)-a(:,id+i)))
    3532                       d(:,jb+j) = s1*((a(:,ia+i)-a(:,ic+i))-(b(:,ib+i)+b(:,id+i))) + c1*((b(:,ia+i)+b(:,ic+i))&
    3533                                 +(a(:,ib+i)-a(:,id+i)))
    3534                       c(:,jd+j) = c3*((a(:,ia+i)-a(:,ic+i))+(b(:,ib+i)+b(:,id+i))) - s3*((b(:,ia+i)+b(:,ic+i))&
    3535                                 -(a(:,ib+i)-a(:,id+i)))
    3536                       d(:,jd+j) = s3*((a(:,ia+i)-a(:,ic+i))+(b(:,ib+i)+b(:,id+i))) + c3*((b(:,ia+i)+b(:,ic+i))&
    3537                                 -(a(:,ib+i)-a(:,id+i)))
     3623                      c(:,ja+j) = ( a(:,ia+i) + a(:,ic+i) ) + ( a(:,ib+i) + a(:,id+i) )
     3624                      d(:,ja+j) = ( b(:,ia+i) - b(:,ic+i) ) + ( b(:,ib+i) - b(:,id+i) )
     3625                      c(:,jc+j) = c2 * ( ( a(:,ia+i) + a(:,ic+i) ) - ( a(:,ib+i) + a(:,id+i) ) )   &
     3626                                  - s2 * ( ( b(:,ia+i) - b(:,ic+i) ) - ( b(:,ib+i) - b(:,id+i) ) )
     3627                      d(:,jc+j) = s2 * ( ( a(:,ia+i) + a(:,ic+i) ) - ( a(:,ib+i) + a(:,id+i) ) )   &
     3628                                  + c2 * ( ( b(:,ia+i) - b(:,ic+i) ) - ( b(:,ib+i) - b(:,id+i) ) )
     3629                      c(:,jb+j) = c1 * ( ( a(:,ia+i) - a(:,ic+i) ) - ( b(:,ib+i) + b(:,id+i) ) )   &
     3630                                  - s1 * ( ( b(:,ia+i) + b(:,ic+i) ) + ( a(:,ib+i) - a(:,id+i) ) )
     3631                      d(:,jb+j) = s1 * ( ( a(:,ia+i) - a(:,ic+i) ) - ( b(:,ib+i) + b(:,id+i) ) )   &
     3632                                  + c1 * ( ( b(:,ia+i) + b(:,ic+i) ) + ( a(:,ib+i) - a(:,id+i) ) )
     3633                      c(:,jd+j) = c3 * ( ( a(:,ia+i) - a(:,ic+i) ) + ( b(:,ib+i) + b(:,id+i) ) )   &
     3634                                  - s3 * ( ( b(:,ia+i) + b(:,ic+i) ) - ( a(:,ib+i) - a(:,id+i) ) )
     3635                      d(:,jd+j) = s3 * ( ( a(:,ia+i) - a(:,ic+i) ) + ( b(:,ib+i) + b(:,id+i) ) )   &
     3636                                  + c3 * ( ( b(:,ia+i) + b(:,ic+i) ) - ( a(:,ib+i) - a(:,id+i) ) )
    35383637                      ibase = ibase + inc1
    35393638                      jbase = jbase + inc2
     
    35563655                j = jbase
    35573656                c(:,ja+j) = a(:,ia+i) + a(:,ib+i)
    3558                 c(:,jb+j) = sin45*((a(:,ia+i)-a(:,ib+i))-(b(:,ia+i)+b(:,ib+i)))
     3657                c(:,jb+j) = sin45 * ( ( a(:,ia+i) - a(:,ib+i) ) - ( b(:,ia+i) + b(:,ib+i) ) )
    35593658                c(:,jc+j) = b(:,ib+i) - b(:,ia+i)
    3560                 c(:,jd+j) = -sin45*((a(:,ia+i)-a(:,ib+i))+(b(:,ia+i)+b(:,ib+i)))
     3659                c(:,jd+j) = - sin45 * ( ( a(:,ia+i) - a(:,ib+i) ) + ( b(:,ia+i) + b(:,ib+i) ) )
    35613660                ibase = ibase + inc1
    35623661                jbase = jbase + inc2
     
    35683667                i = ibase
    35693668                j = jbase
    3570                 c(:,ja+j) = 2.0_wp*((a(:,ia+i)+a(:,ic+i))+a(:,ib+i))
    3571                 c(:,jb+j) = 2.0_wp*((a(:,ia+i)-a(:,ic+i))-b(:,ib+i))
    3572                 c(:,jc+j) = 2.0_wp*((a(:,ia+i)+a(:,ic+i))-a(:,ib+i))
    3573                 c(:,jd+j) = 2.0_wp*((a(:,ia+i)-a(:,ic+i))+b(:,ib+i))
     3669                c(:,ja+j) = 2.0_wp * ( ( a(:,ia+i) + a(:,ic+i) ) + a(:,ib+i) )
     3670                c(:,jb+j) = 2.0_wp * ( ( a(:,ia+i) - a(:,ic+i) ) - b(:,ib+i) )
     3671                c(:,jc+j) = 2.0_wp * ( ( a(:,ia+i) + a(:,ic+i) ) - a(:,ib+i) )
     3672                c(:,jd+j) = 2.0_wp * ( ( a(:,ia+i) - a(:,ic+i) ) + b(:,ib+i) )
    35743673                ibase = ibase + inc1
    35753674                jbase = jbase + inc2
     
    35833682
    35843683          ia = 1
    3585           ib = ia + (2*m-la) * inc1
    3586           ic = ib + 2*m*inc1
     3684          ib = ia + ( 2 * m - la ) * inc1
     3685          ic = ib + 2 * m * inc1
    35873686          id = ic
    35883687          ie = ib
     
    35983697                i = ibase
    35993698                j = jbase
    3600                 c(:,ja+j) = a(:,ia+i) + (a(:,ib+i)+a(:,ic+i))
    3601                 c(:,jb+j) = ((a(:,ia+i)-0.25_wp*(a(:,ib+i)+a(:,ic+i)))+qrt5*(a(:,ib+i)-a(:,ic+i))) -        &
    3602                           (sin72*b(:,ib+i)+sin36*b(:,ic+i))
    3603                 c(:,jc+j) = ((a(:,ia+i)-0.25_wp*(a(:,ib+i)+a(:,ic+i)))-qrt5*(a(:,ib+i)-a(:,ic+i))) -        &
    3604                           (sin36*b(:,ib+i)-sin72*b(:,ic+i))
    3605                 c(:,jd+j) = ((a(:,ia+i)-0.25_wp*(a(:,ib+i)+a(:,ic+i)))-qrt5*(a(:,ib+i)-a(:,ic+i))) +        &
    3606                           (sin36*b(:,ib+i)-sin72*b(:,ic+i))
    3607                 c(:,je+j) = ((a(:,ia+i)-0.25_wp*(a(:,ib+i)+a(:,ic+i)))+qrt5*(a(:,ib+i)-a(:,ic+i))) +        &
    3608                           (sin72*b(:,ib+i)+sin36*b(:,ic+i))
     3699                c(:,ja+j) = a(:,ia+i) + ( a(:,ib+i) + a(:,ic+i) )
     3700                c(:,jb+j) = ( ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) ) + qrt5           &
     3701                          * ( a(:,ib+i) - a(:,ic+i) ) ) - ( sin72 * b(:,ib+i) + sin36 * b(:,ic+i) )
     3702                c(:,jc+j) = ( ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) ) - qrt5           &
     3703                          * ( a(:,ib+i) - a(:,ic+i) ) ) - ( sin36 * b(:,ib+i) - sin72 * b(:,ic+i) )
     3704                c(:,jd+j) = ( ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) ) - qrt5           &
     3705                          * ( a(:,ib+i) - a(:,ic+i) ) ) + ( sin36 * b(:,ib+i) - sin72 * b(:,ic+i) )
     3706                c(:,je+j) = ( ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) ) + qrt5           &
     3707                          * ( a(:,ib+i) - a(:,ic+i) ) ) + ( sin72 * b(:,ib+i) + sin36 * b(:,ic+i) )
    36093708                ibase = ibase + inc1
    36103709                jbase = jbase + inc2
     
    36393738                      j = jbase
    36403739!DIR$ IVDEP
    3641                       DO  ipl = 1, SIZE(a,1)
    3642                          a10      = (a(ipl,ia+i)-0.25_wp*((a(ipl,ib+i)+a(ipl,ie+i))+(a(ipl,ic+i)+a(ipl,id+i)))) +      &
    3643                                     qrt5*((a(ipl,ib+i)+a(ipl,ie+i))-(a(ipl,ic+i)+a(ipl,id+i)))
    3644                          a20      = (a(ipl,ia+i)-0.25_wp*((a(ipl,ib+i)+a(ipl,ie+i))+(a(ipl,ic+i)+a(ipl,id+i)))) -      &
    3645                                     qrt5*((a(ipl,ib+i)+a(ipl,ie+i))-(a(ipl,ic+i)+a(ipl,id+i)))
    3646                          b10      = (b(ipl,ia+i)-0.25_wp*((b(ipl,ib+i)-b(ipl,ie+i))+(b(ipl,ic+i)-b(ipl,id+i)))) +      &
    3647                                     qrt5*((b(ipl,ib+i)-b(ipl,ie+i))-(b(ipl,ic+i)-b(ipl,id+i)))
    3648                          b20      = (b(ipl,ia+i)-0.25_wp*((b(ipl,ib+i)-b(ipl,ie+i))+(b(ipl,ic+i)-b(ipl,id+i)))) -      &
    3649                                     qrt5*((b(ipl,ib+i)-b(ipl,ie+i))-(b(ipl,ic+i)-b(ipl,id+i)))
    3650                          a11      = sin72*(b(ipl,ib+i)+b(ipl,ie+i)) + sin36*(b(ipl,ic+i)+b(ipl,id+i))
    3651                          a21      = sin36*(b(ipl,ib+i)+b(ipl,ie+i)) - sin72*(b(ipl,ic+i)+b(ipl,id+i))
    3652                          b11      = sin72*(a(ipl,ib+i)-a(ipl,ie+i)) + sin36*(a(ipl,ic+i)-a(ipl,id+i))
    3653                          b21      = sin36*(a(ipl,ib+i)-a(ipl,ie+i)) - sin72*(a(ipl,ic+i)-a(ipl,id+i))
    3654                          c(ipl,ja+j) = a(ipl,ia+i) + ((a(ipl,ib+i)+a(ipl,ie+i))+(a(ipl,ic+i)+a(ipl,id+i)))
    3655                          d(ipl,ja+j) = b(ipl,ia+i) + ((b(ipl,ib+i)-b(ipl,ie+i))+(b(ipl,ic+i)-b(ipl,id+i)))
    3656                          c(ipl,jb+j) = c1*(a10     -a11     ) - s1*(b10     +b11     )
    3657                          d(ipl,jb+j) = s1*(a10     -a11     ) + c1*(b10     +b11     )
    3658                          c(ipl,je+j) = c4*(a10     +a11     ) - s4*(b10     -b11     )
    3659                          d(ipl,je+j) = s4*(a10     +a11     ) + c4*(b10     -b11     )
    3660                          c(ipl,jc+j) = c2*(a20     -a21     ) - s2*(b20     +b21     )
    3661                          d(ipl,jc+j) = s2*(a20     -a21     ) + c2*(b20     +b21     )
    3662                          c(ipl,jd+j) = c3*(a20     +a21     ) - s3*(b20     -b21     )
    3663                          d(ipl,jd+j) = s3*(a20     +a21     ) + c3*(b20     -b21     )
     3740                      DO  ipl = 1, SIZE( a, 1 )
     3741                         a10 = ( a (ipl,ia+i) - 0.25_wp * ( ( a(ipl,ib+i) + a(ipl,ie+i) )          &
     3742                             + ( a(ipl,ic+i) + a(ipl,id+i) ) ) ) + qrt5 * ( ( a(ipl,ib+i)          &
     3743                             + a(ipl,ie+i) ) - ( a(ipl,ic+i) + a(ipl,id+i) ) )
     3744                         a20 = ( a(ipl,ia+i) - 0.25_wp * ( ( a(ipl,ib+i) + a(ipl,ie+i) )           &
     3745                             + ( a(ipl,ic+i) + a(ipl,id+i) ) ) ) - qrt5 * ( ( a(ipl,ib+i)          &
     3746                             + a(ipl,ie+i) ) - ( a(ipl,ic+i) + a(ipl,id+i) ) )
     3747                         b10 = ( b(ipl,ia+i) - 0.25_wp * ( ( b(ipl,ib+i) - b(ipl,ie+i) )           &
     3748                             + ( b(ipl,ic+i) - b(ipl,id+i) ) ) ) + qrt5 * ( ( b(ipl,ib+i)          &
     3749                             - b(ipl,ie+i) ) - ( b(ipl,ic+i) - b(ipl,id+i) ) )
     3750                         b20 = ( b(ipl,ia+i) - 0.25_wp * ( ( b(ipl,ib+i) - b(ipl,ie+i) )           &
     3751                             + ( b(ipl,ic+i) - b(ipl,id+i) ) ) ) - qrt5 * ( ( b(ipl,ib+i)          &
     3752                             - b(ipl,ie+i) ) - ( b(ipl,ic+i) - b(ipl,id+i) ) )
     3753                         a11 = sin72 * ( b(ipl,ib+i) + b(ipl,ie+i) ) + sin36                       &
     3754                               * ( b(ipl,ic+i) + b(ipl,id+i) )
     3755                         a21 = sin36 * ( b(ipl,ib+i) + b(ipl,ie+i) ) - sin72                       &
     3756                               * ( b(ipl,ic+i) + b(ipl,id+i) )
     3757                         b11 = sin72 * ( a(ipl,ib+i) - a(ipl,ie+i) ) + sin36                       &
     3758                               * ( a(ipl,ic+i) - a(ipl,id+i) )
     3759                         b21 = sin36 * ( a(ipl,ib+i) - a(ipl,ie+i) ) - sin72                       &
     3760                               * ( a(ipl,ic+i) - a(ipl,id+i) )
     3761                         c(ipl,ja+j) = a(ipl,ia+i) + ( ( a(ipl,ib+i) + a(ipl,ie+i) )               &
     3762                                       + ( a(ipl,ic+i) + a(ipl,id+i) ) )
     3763                         d(ipl,ja+j) = b(ipl,ia+i) + ( ( b(ipl,ib+i) - b(ipl,ie+i) )               &
     3764                                       + ( b(ipl,ic+i) - b(ipl,id+i) ) )
     3765                         c(ipl,jb+j) = c1 * ( a10 - a11 ) - s1 * ( b10 + b11 )
     3766                         d(ipl,jb+j) = s1 * ( a10 - a11 ) + c1 * ( b10 + b11 )
     3767                         c(ipl,je+j) = c4 * ( a10 + a11 ) - s4 * ( b10 - b11 )
     3768                         d(ipl,je+j) = s4 * ( a10 + a11 ) + c4 * ( b10 - b11 )
     3769                         c(ipl,jc+j) = c2 * ( a20 - a21 ) - s2 * ( b20 + b21 )
     3770                         d(ipl,jc+j) = s2 * ( a20 - a21 ) + c2 * ( b20 + b21 )
     3771                         c(ipl,jd+j) = c3 * ( a20 + a21 ) - s3 * ( b20 - b21 )
     3772                         d(ipl,jd+j) = s3 * ( a20 + a21 ) + c3 * ( b20 - b21 )
    36643773                      ENDDO
    36653774                      ibase = ibase + inc1
     
    36823791                i = ibase
    36833792                j = jbase
    3684                 c(:,ja+j) = (a(:,ia+i)+a(:,ib+i)) + a(:,ic+i)
    3685                 c(:,jb+j) = (qrt5*(a(:,ia+i)-a(:,ib+i))+(0.25_wp*(a(:,ia+i)+a(:,ib+i))-a(:,ic+i))) -        &
    3686                           (sin36*b(:,ia+i)+sin72*b(:,ib+i))
    3687                 c(:,je+j) = -(qrt5*(a(:,ia+i)-a(:,ib+i))+(0.25_wp*(a(:,ia+i)+a(:,ib+i))-a(:,ic+i))) -       &
    3688                           (sin36*b(:,ia+i)+sin72*b(:,ib+i))
    3689                 c(:,jc+j) = (qrt5*(a(:,ia+i)-a(:,ib+i))-(0.25_wp*(a(:,ia+i)+a(:,ib+i))-a(:,ic+i))) -        &
    3690                           (sin72*b(:,ia+i)-sin36*b(:,ib+i))
    3691                 c(:,jd+j) = -(qrt5*(a(:,ia+i)-a(:,ib+i))-(0.25_wp*(a(:,ia+i)+a(:,ib+i))-a(:,ic+i))) -       &
    3692                           (sin72*b(:,ia+i)-sin36*b(:,ib+i))
     3793                c(:,ja+j) = ( a(:,ia+i) + a(:,ib+i) ) + a(:,ic+i)
     3794                c(:,jb+j) = ( qrt5 * ( a(:,ia+i) - a(:,ib+i) ) + ( 0.25_wp * ( a(:,ia+i)           &
     3795                            + a(:,ib+i) ) - a(:,ic+i) ) ) - ( sin36 * b(:,ia+i) + sin72 * b(:,ib+i) )
     3796                c(:,je+j) = - ( qrt5 * ( a(:,ia+i) - a(:,ib+i) ) + ( 0.25_wp * ( a(:,ia+i)         &
     3797                            + a(:,ib+i) ) - a(:,ic+i) ) ) - ( sin36 * b(:,ia+i) + sin72 * b(:,ib+i) )
     3798                c(:,jc+j) = ( qrt5 * ( a(:,ia+i) - a(:,ib+i) ) - ( 0.25_wp * ( a(:,ia+i)           &
     3799                            + a(:,ib+i) ) - a(:,ic+i) ) ) - ( sin72 * b(:,ia+i) - sin36 * b(:,ib+i) )
     3800                c(:,jd+j) = - ( qrt5 * ( a(:,ia+i) - a(:,ib+i) ) - ( 0.25_wp * ( a(:,ia+i)         &
     3801                            + a(:,ib+i) ) - a(:,ic+i) ) ) - ( sin72 * b(:,ia+i) - sin36 * b(:,ib+i) )
    36933802                ibase = ibase + inc1
    36943803                jbase = jbase + inc2
     
    37033812                i = ibase
    37043813                j = jbase
    3705                 c(:,ja+j) = 2.0_wp*(a(:,ia+i)+(a(:,ib+i)+a(:,ic+i)))
    3706                 c(:,jb+j) = (2.0_wp*(a(:,ia+i)-0.25_wp*(a(:,ib+i)+a(:,ic+i)))+qqrt5*(a(:,ib+i)-a(:,ic+i)))  &
    3707                           - (ssin72*b(:,ib+i)+ssin36*b(:,ic+i))
    3708                 c(:,jc+j) = (2.0_wp*(a(:,ia+i)-0.25_wp*(a(:,ib+i)+a(:,ic+i)))-qqrt5*(a(:,ib+i)-a(:,ic+i)))  &
    3709                           - (ssin36*b(:,ib+i)-ssin72*b(:,ic+i))
    3710                 c(:,jd+j) = (2.0_wp*(a(:,ia+i)-0.25_wp*(a(:,ib+i)+a(:,ic+i)))-qqrt5*(a(:,ib+i)-a(:,ic+i)))  &
    3711                           + (ssin36*b(:,ib+i)-ssin72*b(:,ic+i))
    3712                 c(:,je+j) = (2.0_wp*(a(:,ia+i)-0.25_wp*(a(:,ib+i)+a(:,ic+i)))+qqrt5*(a(:,ib+i)-a(:,ic+i)))  &
    3713                           + (ssin72*b(:,ib+i)+ssin36*b(:,ic+i))
     3814                c(:,ja+j) = 2.0_wp * ( a(:,ia+i) + ( a(:,ib+i) + a(:,ic+i) ) )
     3815                c(:,jb+j) = ( 2.0_wp * ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) )         &
     3816                            + qqrt5 * ( a(:,ib+i) - a(:,ic+i) ) )  - ( ssin72 * b(:,ib+i)          &
     3817                            + ssin36 * b(:,ic+i) )
     3818                c(:,jc+j) = ( 2.0_wp * ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) )         &
     3819                            - qqrt5 * ( a(:,ib+i) - a(:,ic+i) ) ) - ( ssin36 * b(:,ib+i)           &
     3820                            - ssin72 * b(:,ic+i) )
     3821                c(:,jd+j) = ( 2.0_wp * ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) )         &
     3822                            - qqrt5 * ( a(:,ib+i) - a(:,ic+i) ) ) + ( ssin36 * b(:,ib+i)           &
     3823                            - ssin72 * b(:,ic+i) )
     3824                c(:,je+j) = ( 2.0_wp * ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) )         &
     3825                            + qqrt5 * ( a(:,ib+i) - a(:,ic+i) ) ) + ( ssin72 * b(:,ib+i)           &
     3826                            + ssin36 * b(:,ic+i) )
    37143827                ibase = ibase + inc1
    37153828                jbase = jbase + inc2
     
    37233836
    37243837          ia = 1
    3725           ib = ia + (2*m-la) * inc1
    3726           ic = ib + 2*m*inc1
    3727           id = ic + 2*m*inc1
     3838          ib = ia + ( 2 * m - la ) * inc1
     3839          ic = ib + 2 * m * inc1
     3840          id = ic + 2 * m * inc1
    37283841          ie = ic
    37293842          if = ib
     
    37403853                i = ibase
    37413854                j = jbase
    3742                 c(:,ja+j) = (a(:,ia+i)+a(:,id+i)) + (a(:,ib+i)+a(:,ic+i))
    3743                 c(:,jd+j) = (a(:,ia+i)-a(:,id+i)) - (a(:,ib+i)-a(:,ic+i))
    3744                 c(:,jb+j) = ((a(:,ia+i)-a(:,id+i))+0.5_wp*(a(:,ib+i)-a(:,ic+i))) - (sin60*(b(:,ib+i)+b(:,ic+i)))
    3745                 c(:,jf+j) = ((a(:,ia+i)-a(:,id+i))+0.5_wp*(a(:,ib+i)-a(:,ic+i))) + (sin60*(b(:,ib+i)+b(:,ic+i)))
    3746                 c(:,jc+j) = ((a(:,ia+i)+a(:,id+i))-0.5_wp*(a(:,ib+i)+a(:,ic+i))) - (sin60*(b(:,ib+i)-b(:,ic+i)))
    3747                 c(:,je+j) = ((a(:,ia+i)+a(:,id+i))-0.5_wp*(a(:,ib+i)+a(:,ic+i))) + (sin60*(b(:,ib+i)-b(:,ic+i)))
     3855                c(:,ja+j) = ( a(:,ia+i) + a(:,id+i) ) +  ( a(:,ib+i) + a(:,ic+i) )
     3856                c(:,jd+j) = ( a(:,ia+i) - a(:,id+i) ) -  ( a(:,ib+i) - a(:,ic+i) )
     3857                c(:,jb+j) = ( ( a(:,ia+i) - a(:,id+i) ) + 0.5_wp * ( a(:,ib+i) - a(:,ic+i) ) )     &
     3858                          - ( sin60 * ( b(:,ib+i) + b(:,ic+i) ) )
     3859                c(:,jf+j) = ( ( a(:,ia+i) - a(:,id+i) ) + 0.5_wp * ( a(:,ib+i) - a(:,ic+i) ) )     &
     3860                          + ( sin60 * ( b(:,ib+i) + b(:,ic+i) ) )
     3861                c(:,jc+j) = ( ( a(:,ia+i) + a(:,id+i) ) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) )     &
     3862                          - ( sin60 * ( b(:,ib+i) - b(:,ic+i) ) )
     3863                c(:,je+j) = ( ( a(:,ia+i) + a(:,id+i) ) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) )     &
     3864                          + ( sin60 * ( b(:,ib+i) - b(:,ic+i) ) )
    37483865                ibase = ibase + inc1
    37493866                jbase = jbase + inc2
     
    37813898                      i = ibase
    37823899                      j = jbase
    3783                       DO  ipl = 1, SIZE(a,1)
    3784                          a11      = (a(ipl,ie+i)+a(ipl,ib+i)) + (a(ipl,ic+i)+a(ipl,if+i))
    3785                          a20      = (a(ipl,ia+i)+a(ipl,id+i)) - 0.5_wp*a11
    3786                          a21      = sin60*((a(ipl,ie+i)+a(ipl,ib+i))-(a(ipl,ic+i)+a(ipl,if+i)))
    3787                          b11      = (b(ipl,ib+i)-b(ipl,ie+i)) + (b(ipl,ic+i)-b(ipl,if+i))
    3788                          b20      = (b(ipl,ia+i)-b(ipl,id+i)) - 0.5_wp*b11
    3789                          b21      = sin60*((b(ipl,ib+i)-b(ipl,ie+i))-(b(ipl,ic+i)-b(ipl,if+i)))
    3790 
    3791                          c(ipl,ja+j) = (a(ipl,ia+i)+a(ipl,id+i)) + a11
    3792                          d(ipl,ja+j) = (b(ipl,ia+i)-b(ipl,id+i)) + b11
    3793                          c(ipl,jc+j) = c2*(a20     -b21     ) - s2*(b20     +a21     )
    3794                          d(ipl,jc+j) = s2*(a20     -b21     ) + c2*(b20     +a21     )
    3795                          c(ipl,je+j) = c4*(a20     +b21     ) - s4*(b20     -a21     )
    3796                          d(ipl,je+j) = s4*(a20     +b21     ) + c4*(b20     -a21     )
    3797 
    3798                          a11      = (a(ipl,ie+i)-a(ipl,ib+i)) + (a(ipl,ic+i)-a(ipl,if+i))
    3799                          b11      = (b(ipl,ie+i)+b(ipl,ib+i)) - (b(ipl,ic+i)+b(ipl,if+i))
    3800                          a20      = (a(ipl,ia+i)-a(ipl,id+i)) - 0.5_wp*a11
    3801                          a21      = sin60*((a(ipl,ie+i)-a(ipl,ib+i))-(a(ipl,ic+i)-a(ipl,if+i)))
    3802                          b20      = (b(ipl,ia+i)+b(ipl,id+i)) + 0.5_wp*b11
    3803                          b21      = sin60*((b(ipl,ie+i)+b(ipl,ib+i))+(b(ipl,ic+i)+b(ipl,if+i)))
    3804 
    3805                          c(ipl,jd+j) = c3*((a(ipl,ia+i)-a(ipl,id+i))+a11     ) - s3*((b(ipl,ia+i)+b(ipl,id+i))-b11     )
    3806                          d(ipl,jd+j) = s3*((a(ipl,ia+i)-a(ipl,id+i))+a11     ) + c3*((b(ipl,ia+i)+b(ipl,id+i))-b11     )
    3807                          c(ipl,jb+j) = c1*(a20     -b21     ) - s1*(b20     -a21     )
    3808                          d(ipl,jb+j) = s1*(a20     -b21     ) + c1*(b20     -a21     )
    3809                          c(ipl,jf+j) = c5*(a20     +b21     ) - s5*(b20     +a21     )
    3810                          d(ipl,jf+j) = s5*(a20     +b21     ) + c5*(b20     +a21     )
     3900                      DO  ipl = 1, SIZE( a, 1)
     3901                         a11 = ( a(ipl,ie+i) + a(ipl,ib+i) ) + ( a(ipl,ic+i) + a(ipl,if+i) )
     3902                         a20 = ( a(ipl,ia+i) + a(ipl,id+i) ) - 0.5_wp * a11
     3903                         a21 = sin60 * ( ( a(ipl,ie+i) + a(ipl,ib+i) )                             &
     3904                               - ( a(ipl,ic+i) + a(ipl,if+i) ) )
     3905                         b11 = ( b(ipl,ib+i) - b(ipl,ie+i) ) + ( b(ipl,ic+i) - b(ipl,if+i) )
     3906                         b20 = ( b(ipl,ia+i) - b(ipl,id+i) ) - 0.5_wp * b11
     3907                         b21 = sin60 * ( ( b(ipl,ib+i) - b(ipl,ie+i) )                             &
     3908                               - ( b(ipl,ic+i) - b(ipl,if+i) ) )
     3909
     3910                         c(ipl,ja+j) = ( a(ipl,ia+i) + a(ipl,id+i) ) + a11
     3911                         d(ipl,ja+j) = ( b(ipl,ia+i) - b(ipl,id+i) ) + b11
     3912                         c(ipl,jc+j) = c2 * ( a20 - b21 ) - s2 * ( b20 + a21 )
     3913                         d(ipl,jc+j) = s2 * ( a20 - b21 ) + c2 * ( b20 + a21 )
     3914                         c(ipl,je+j) = c4 * ( a20 + b21 ) - s4 * ( b20 - a21 )
     3915                         d(ipl,je+j) = s4 * ( a20 + b21 ) + c4 * ( b20 - a21 )
     3916
     3917                         a11 = ( a(ipl,ie+i) - a(ipl,ib+i) ) + ( a(ipl,ic+i) - a(ipl,if+i) )
     3918                         b11 = ( b(ipl,ie+i) + b(ipl,ib+i) ) - ( b(ipl,ic+i) + b(ipl,if+i) )
     3919                         a20 = ( a(ipl,ia+i) - a(ipl,id+i) ) - 0.5_wp * a11
     3920                         a21 = sin60 * ( ( a(ipl,ie+i) - a(ipl,ib+i) ) - ( a(ipl,ic+i)             &
     3921                               - a(ipl,if+i) ) )
     3922                         b20 = ( b(ipl,ia+i) + b(ipl,id+i) ) + 0.5_wp * b11
     3923                         b21 = sin60 * ( ( b(ipl,ie+i) + b(ipl,ib+i) ) + ( b(ipl,ic+i)             &
     3924                               + b(ipl,if+i) ) )
     3925
     3926                         c(ipl,jd+j) = c3 * ( ( a(ipl,ia+i) - a(ipl,id+i) ) + a11 )                &
     3927                                     - s3 * ( ( b(ipl,ia+i) + b(ipl,id+i) ) - b11 )
     3928                         d(ipl,jd+j) = s3 * ( ( a(ipl,ia+i) - a(ipl,id+i) ) + a11 )                &
     3929                                     + c3 * ( ( b(ipl,ia+i) + b(ipl,id+i) ) - b11 )
     3930                         c(ipl,jb+j) = c1 * ( a20 - b21 ) - s1 * ( b20 - a21 )
     3931                         d(ipl,jb+j) = s1 * ( a20 - b21 ) + c1 * ( b20 - a21 )
     3932                         c(ipl,jf+j) = c5 * ( a20 + b21 ) - s5 * ( b20 + a21 )
     3933                         d(ipl,jf+j) = s5 * ( a20 + b21 ) + c5 * ( b20 + a21 )
    38113934
    38123935                      ENDDO
     
    38313954                i = ibase
    38323955                j = jbase
    3833                 DO  ipl = 1, SIZE(a,1)
    3834                    c(ipl,ja+j) = a(ipl,ib+i) + (a(ipl,ia+i)+a(ipl,ic+i))
    3835                    c(ipl,jd+j) = b(ipl,ib+i) - (b(ipl,ia+i)+b(ipl,ic+i))
    3836                    c(ipl,jb+j) = (sin60*(a(ipl,ia+i)-a(ipl,ic+i))) - (0.5_wp*(b(ipl,ia+i)+b(ipl,ic+i))+b(ipl,ib+i))
    3837                    c(ipl,jf+j) = -(sin60*(a(ipl,ia+i)-a(ipl,ic+i))) - (0.5_wp*(b(ipl,ia+i)+b(ipl,ic+i))+b(ipl,ib+i))
    3838                    c(ipl,jc+j) = sin60*(b(ipl,ic+i)-b(ipl,ia+i)) + (0.5_wp*(a(ipl,ia+i)+a(ipl,ic+i))-a(ipl,ib+i))
    3839                    c(ipl,je+j) = sin60*(b(ipl,ic+i)-b(ipl,ia+i)) - (0.5_wp*(a(ipl,ia+i)+a(ipl,ic+i))-a(ipl,ib+i))
     3956                DO  ipl = 1, SIZE( a ,1 )
     3957                   c(ipl,ja+j) = a(ipl,ib+i) + ( a(ipl,ia+i) + a(ipl,ic+i) )
     3958                   c(ipl,jd+j) = b(ipl,ib+i) - ( b(ipl,ia+i) + b(ipl,ic+i) )
     3959                   c(ipl,jb+j) = ( sin60 * ( a(ipl,ia+i) - a(ipl,ic+i) ) )                         &
     3960                               - ( 0.5_wp * ( b(ipl,ia+i) + b(ipl,ic+i) ) + b(ipl,ib+i) )
     3961                   c(ipl,jf+j) = - ( sin60 * ( a(ipl,ia+i) - a(ipl,ic+i) ) )                       &
     3962                                 - ( 0.5_wp * ( b(ipl,ia+i) + b(ipl,ic+i) ) + b(ipl,ib+i) )
     3963                   c(ipl,jc+j) = sin60 * ( b(ipl,ic+i) - b(ipl,ia+i) )                              &
     3964                                 + ( 0.5_wp * ( a(ipl,ia+i) + a(ipl,ic+i) ) - a(ipl,ib+i) )
     3965                   c(ipl,je+j) = sin60 * ( b(ipl,ic+i) - b(ipl,ia+i) )                              &
     3966                                 - ( 0.5_wp * ( a(ipl,ia+i) + a(ipl,ic+i) ) - a(ipl,ib+i) )
    38403967                ENDDO
    38413968                ibase = ibase + inc1
     
    38493976                i = ibase
    38503977                j = jbase
    3851                 c(:,ja+j) = (2.0_wp*(a(:,ia+i)+a(:,id+i))) + (2.0_wp*(a(:,ib+i)+a(:,ic+i)))
    3852                 c(:,jd+j) = (2.0_wp*(a(:,ia+i)-a(:,id+i))) - (2.0_wp*(a(:,ib+i)-a(:,ic+i)))
    3853                 c(:,jb+j) = (2.0_wp*(a(:,ia+i)-a(:,id+i))+(a(:,ib+i)-a(:,ic+i))) - (ssin60*(b(:,ib+i)+b(:,ic+i)))
    3854                 c(:,jf+j) = (2.0_wp*(a(:,ia+i)-a(:,id+i))+(a(:,ib+i)-a(:,ic+i))) + (ssin60*(b(:,ib+i)+b(:,ic+i)))
    3855                 c(:,jc+j) = (2.0_wp*(a(:,ia+i)+a(:,id+i))-(a(:,ib+i)+a(:,ic+i))) - (ssin60*(b(:,ib+i)-b(:,ic+i)))
    3856                 c(:,je+j) = (2.0_wp*(a(:,ia+i)+a(:,id+i))-(a(:,ib+i)+a(:,ic+i))) + (ssin60*(b(:,ib+i)-b(:,ic+i)))
     3978                c(:,ja+j) = ( 2.0_wp * ( a(:,ia+i) + a(:,id+i) ) )                                 &
     3979                          + ( 2.0_wp * ( a(:,ib+i) + a(:,ic+i) ) )
     3980                c(:,jd+j) = ( 2.0_wp * ( a(:,ia+i) - a(:,id+i) ) )                                 &
     3981                          - ( 2.0_wp * ( a(:,ib+i) - a(:,ic+i) ) )
     3982                c(:,jb+j) = ( 2.0_wp * ( a(:,ia+i) - a(:,id+i) ) + ( a(:,ib+i) - a(:,ic+i) ) )     &
     3983                          - ( ssin60 * ( b(:,ib+i) + b(:,ic+i) ) )
     3984                c(:,jf+j) = ( 2.0_wp * ( a(:,ia+i) - a(:,id+i) ) + ( a(:,ib+i) - a(:,ic+i) ) )     &
     3985                          + ( ssin60 * ( b(:,ib+i) + b(:,ic+i) ) )
     3986                c(:,jc+j) = ( 2.0_wp * ( a(:,ia+i) + a(:,id+i) ) - ( a(:,ib+i) + a(:,ic+i) ) )     &
     3987                          - ( ssin60 * ( b(:,ib+i) - b(:,ic+i) ) )
     3988                c(:,je+j) = ( 2.0_wp * ( a(:,ia+i) + a(:,id+i) ) - ( a(:,ib+i) + a(:,ic+i) ) )     &
     3989                          + ( ssin60 * ( b(:,ib+i) - b(:,ic+i) ) )
    38573990                ibase = ibase + inc1
    38583991                jbase = jbase + inc2
     
    38714004
    38724005          ia = 1
    3873           ib = ia + la*inc1
    3874           ic = ib + 2*la*inc1
    3875           id = ic + 2*la*inc1
    3876           ie = id + 2*la*inc1
     4006          ib = ia + la * inc1
     4007          ic = ib + 2 * la * inc1
     4008          id = ic + 2 * la * inc1
     4009          ie = id + 2 * la * inc1
    38774010          ja = 1
    38784011          jb = ja + jink
     
    38884021             i = ibase
    38894022             j = jbase
    3890              c(:,ja+j) = 2.0_wp*(((a(:,ia+i)+a(:,ie+i))+a(:,ic+i))+(a(:,ib+i)+a(:,id+i)))
    3891              c(:,je+j) = 2.0_wp*(((a(:,ia+i)+a(:,ie+i))+a(:,ic+i))-(a(:,ib+i)+a(:,id+i)))
    3892              c(:,jc+j) = 2.0_wp*(((a(:,ia+i)+a(:,ie+i))-a(:,ic+i))-(b(:,ib+i)-b(:,id+i)))
    3893              c(:,jg+j) = 2.0_wp*(((a(:,ia+i)+a(:,ie+i))-a(:,ic+i))+(b(:,ib+i)-b(:,id+i)))
    3894              c(:,jb+j) = 2.0_wp*((a(:,ia+i)-a(:,ie+i))-b(:,ic+i)) + ssin45*((a(:,ib+i)-a(:,id+i))-(b(:,ib+i)+b(:,id+i)))
    3895              c(:,jf+j) = 2.0_wp*((a(:,ia+i)-a(:,ie+i))-b(:,ic+i)) - ssin45*((a(:,ib+i)-a(:,id+i))-(b(:,ib+i)+b(:,id+i)))
    3896              c(:,jd+j) = 2.0_wp*((a(:,ia+i)-a(:,ie+i))+b(:,ic+i)) - ssin45*((a(:,ib+i)-a(:,id+i))+(b(:,ib+i)+b(:,id+i)))
    3897              c(:,jh+j) = 2.0_wp*((a(:,ia+i)-a(:,ie+i))+b(:,ic+i)) + ssin45*((a(:,ib+i)-a(:,id+i))+(b(:,ib+i)+b(:,id+i)))
     4023             c(:,ja+j) = 2.0_wp * ( ( ( a(:,ia+i) + a(:,ie+i) ) + a(:,ic+i) )                      &
     4024                                                  + ( a(:,ib+i) + a(:,id+i) ) )
     4025             c(:,je+j) = 2.0_wp * ( ( ( a(:,ia+i) + a(:,ie+i) ) + a(:,ic+i) )                      &
     4026                                                - ( a(:,ib+i) + a(:,id+i) ) )
     4027             c(:,jc+j) = 2.0_wp * ( ( ( a(:,ia+i) + a(:,ie+i) ) - a(:,ic+i) )                      &
     4028                                                  - ( b(:,ib+i) - b(:,id+i) ) )
     4029             c(:,jg+j) = 2.0_wp * ( ( ( a(:,ia+i) + a(:,ie+i) ) - a(:,ic+i) )                      &
     4030                                                  + ( b(:,ib+i) - b(:,id+i) ) )
     4031             c(:,jb+j) = 2.0_wp * ( ( a(:,ia+i) - a(:,ie+i) ) - b(:,ic+i) ) + ssin45               &
     4032                                * ( (a(:,ib+i) - a(:,id+i) ) - ( b(:,ib+i) + b(:,id+i) ) )
     4033             c(:,jf+j) = 2.0_wp * ( ( a(:,ia+i) - a(:,ie+i) ) - b(:,ic+i) ) - ssin45               &
     4034                                * ( (a(:,ib+i) - a(:,id+i) ) - ( b(:,ib+i) + b(:,id+i) ) )
     4035             c(:,jd+j) = 2.0_wp * ( ( a(:,ia+i) - a(:,ie+i) ) + b(:,ic+i) ) - ssin45               &
     4036                                * ( (a(:,ib+i) - a(:,id+i) ) + ( b(:,ib+i) + b(:,id+i) ) )
     4037             c(:,jh+j) = 2.0_wp * ( ( a(:,ia+i) - a(:,ie+i) ) + b(:,ic+i) ) + ssin45               &
     4038                                * ( (a(:,ib+i) - a(:,id+i) ) + ( b(:,ib+i) + b(:,id+i) ) )
    38984039             ibase = ibase + inc1
    38994040             jbase = jbase + inc2
Note: See TracChangeset for help on using the changeset viewer.