Index: palm/trunk/SOURCE/temperton_fft_mod.f90
===================================================================
--- palm/trunk/SOURCE/temperton_fft_mod.f90	(revision 4370)
+++ palm/trunk/SOURCE/temperton_fft_mod.f90	(revision 4559)
@@ -1,4 +1,20 @@
 !> @file temperton_fft_mod.f90
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
+! This file is part of the PALM model system.
+!
+! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
+! Public License as published by the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
+! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
+! Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along with PALM. If not, see
+! <http://www.gnu.org/licenses/>.
+!
+! Copyright 1997-2020 Leibniz Universitaet Hannover
+!--------------------------------------------------------------------------------------------------!
+!
 !
 ! Current revisions:
@@ -9,12 +25,15 @@
 ! -----------------
 ! $Id$
+! File re-formatted to follow the PALM coding standard
+!
+! 4370 2020-01-10 14:00:44Z raasch
 ! unused variables removed
-! 
+!
 ! 4366 2020-01-09 08:12:43Z raasch
 ! vectorized routines added
-! 
+!
 ! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
-! 
+!
 ! 3725 2019-02-07 10:11:02Z raasch
 ! GOTO statements replaced, file re-formatted corresponding to coding standards
@@ -27,5 +46,5 @@
 ! ------------
 !> Fast Fourier transformation developed by Clive Temperton, ECMWF.
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  MODULE temperton_fft
 
@@ -57,17 +76,17 @@
 CONTAINS
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Calls Fortran-versions of fft's.
-!> 
+!>
 !> Method:
-!> 
-!> Subroutine 'fft991cy' - multiple fast real periodic transform
-!> supersedes previous routine 'fft991cy'.
-!> 
-!> Real transform of length n performed by removing redundant
-!> operations from complex transform of length n.
-!> 
+!>
+!> Subroutine 'fft991cy' - multiple fast real periodic transform supersedes previous routine
+!> 'fft991cy'.
+!>
+!> Real transform of length n performed by removing redundant operations from complex transform of
+!> length n.
+!>
 !> a       is the array containing input & output data.
 !> work    is an area of size (n+1)*min(lot,nfft).
@@ -81,28 +100,27 @@
 !> isign = +1 for transform from spectral to gridpoint
 !>       = -1 for transform from gridpoint to spectral.
-!> 
+!>
 !> ordering of coefficients:
 !> a(0),b(0),a(1),b(1),a(2),b(2),.,a(n/2),b(n/2)
 !> where b(0)=b(n/2)=0; (n+2) locations required.
-!> 
+!>
 !> ordering of data:
 !> x(0),x(1),x(2),.,x(n-1), 0 , 0 ; (n+2) locations required.
-!> 
-!> Vectorization is achieved on cray by doing the transforms
-!> in parallel.
-!> 
+!>
+!> Vectorization is achieved on cray by doing the transforms in parallel.
+!>
 !> n must be composed of factors 2,3 & 5 but does not have to be even.
-!> 
+!>
 !> definition of transforms:
-!> 
+!>
 !> isign=+1: x(j)=sum(k=0,.,n-1)(c(k)*exp(2*i*j*k*pi/n))
 !> where c(k)=a(k)+i*b(k) and c(n-k)=a(k)-i*b(k)
-!> 
+!>
 !> isign=-1: a(k)=(1/n)*sum(j=0,.,n-1)(x(j)*cos(2*j*k*pi/n))
 !> b(k)=-(1/n)*sum(j=0,.,n-1)(x(j)*sin(2*j*k*pi/n))
-!> 
+!>
 !> calls fortran-versions of fft's  !!!
 !> dimension a(n),work(n),trigs(n),ifax(1)
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE fft991cy( a, work, trigs, ifax, inc, jump, n, lot, isign )
 
@@ -111,35 +129,35 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  inc   !< 
-    INTEGER(iwp) ::  isign !< 
-    INTEGER(iwp) ::  jump  !< 
-    INTEGER(iwp) ::  lot   !< 
-    INTEGER(iwp) ::  n     !< 
-
-    REAL(wp)     ::  a(*)     !< 
-    REAL(wp)     ::  trigs(*) !< 
-    REAL(wp)     ::  work(*)  !< 
-    INTEGER(iwp) ::  ifax(*)  !< 
-
-    INTEGER(iwp) ::  i      !< 
-    INTEGER(iwp) ::  ia     !< 
-    INTEGER(iwp) ::  ibase  !< 
-    INTEGER(iwp) ::  ierr   !< 
-    INTEGER(iwp) ::  ifac   !< 
-    INTEGER(iwp) ::  igo    !< 
-    INTEGER(iwp) ::  ii     !< 
-    INTEGER(iwp) ::  istart !< 
-    INTEGER(iwp) ::  ix     !< 
-    INTEGER(iwp) ::  iz     !< 
-    INTEGER(iwp) ::  j      !< 
-    INTEGER(iwp) ::  jbase  !< 
-    INTEGER(iwp) ::  jj     !< 
-    INTEGER(iwp) ::  k      !< 
-    INTEGER(iwp) ::  la     !< 
-    INTEGER(iwp) ::  nb     !< 
-    INTEGER(iwp) ::  nblox  !< 
-    INTEGER(iwp) ::  nfax   !< 
-    INTEGER(iwp) ::  nvex   !< 
-    INTEGER(iwp) ::  nx     !< 
+    INTEGER(iwp) ::  inc    !<
+    INTEGER(iwp) ::  isign  !<
+    INTEGER(iwp) ::  jump   !<
+    INTEGER(iwp) ::  lot    !<
+    INTEGER(iwp) ::  n      !<
+
+    INTEGER(iwp) ::  i        !<
+    INTEGER(iwp) ::  ia       !<
+    INTEGER(iwp) ::  ibase    !<
+    INTEGER(iwp) ::  ierr     !<
+    INTEGER(iwp) ::  ifac     !<
+    INTEGER(iwp) ::  igo      !<
+    INTEGER(iwp) ::  ii       !<
+    INTEGER(iwp) ::  ifax(*)  !<
+    INTEGER(iwp) ::  istart   !<
+    INTEGER(iwp) ::  ix       !<
+    INTEGER(iwp) ::  iz       !<
+    INTEGER(iwp) ::  j        !<
+    INTEGER(iwp) ::  jbase    !<
+    INTEGER(iwp) ::  jj       !<
+    INTEGER(iwp) ::  k        !<
+    INTEGER(iwp) ::  la       !<
+    INTEGER(iwp) ::  nb       !<
+    INTEGER(iwp) ::  nblox    !<
+    INTEGER(iwp) ::  nfax     !<
+    INTEGER(iwp) ::  nvex     !<
+    INTEGER(iwp) ::  nx       !<
+
+    REAL(wp)     ::  a(*)      !<
+    REAL(wp)     ::  trigs(*)  !<
+    REAL(wp)     ::  work(*)   !<
 
 
@@ -147,7 +165,7 @@
     nfax = ifax(1)
     nx   = n + 1
-    IF ( MOD(n,2) == 1 )  nx = n
+    IF ( MOD( n, 2) == 1 )  nx = n
     nblox = 1 + (lot-1) / nfft
-    nvex = lot - (nblox-1) * nfft
+    nvex = lot - ( nblox - 1 ) * nfft
 
 
@@ -160,5 +178,5 @@
 
           ia = istart
-          i = istart
+          i  = istart
           !DIR$ IVDEP
           !CDIR NODEP
@@ -166,11 +184,11 @@
           DO  j = 1, nvex
              a(i+inc) = 0.5_wp * a(i)
-             i = i + jump
+             i        = i + jump
           ENDDO
           IF ( MOD(n,2) /= 1 )  THEN
              i = istart + n * inc
-             DO  j = 1, nvex
+             DO  j   = 1, nvex
                 a(i) = 0.5_wp * a(i)
-                i = i + jump
+                i    = i + jump
              ENDDO
           ENDIF
@@ -211,5 +229,5 @@
 !
 !--       If necessary, copy results back to a
-          IF ( MOD(nfax,2) /= 0 )  THEN
+          IF ( MOD( nfax, 2 ) /= 0 )  THEN
              ibase = 1
              jbase = ia
@@ -219,6 +237,6 @@
                 DO  ii = 1, n
                    a(j) = work(i)
-                   i = i + 1
-                   j = j + inc
+                   i    = i + 1
+                   j    = j + inc
                 ENDDO
                 ibase = ibase + nx
@@ -228,14 +246,14 @@
 !
 !--       Fill in zeros at end
-          ix = istart + n*inc
+          ix = istart + n * inc
           !DIR$ IVDEP
           !CDIR NODEP
           !OCL NOVREC
           DO  j = 1, nvex
-             a(ix) = 0.0_wp
+             a(ix)     = 0.0_wp
              a(ix+inc) = 0.0_wp
-             ix = ix + jump
+             ix        = ix + jump
           ENDDO
-          istart = istart + nvex*jump
+          istart = istart + nvex * jump
           nvex = nfft
 
@@ -254,5 +272,5 @@
 
              ifac = ifax(nfax+2-k)
-             la = la / ifac
+             la   = la / ifac
              ierr = -1
              IF ( igo /= -1 )  THEN
@@ -291,6 +309,6 @@
                 DO  ii = 1, n
                    a(j) = work(i)
-                   i = i + 1
-                   j = j + inc
+                   i    = i + 1
+                   j    = j + inc
                 ENDDO
                 ibase = ibase + nx
@@ -305,7 +323,7 @@
           !OCL NOVREC
           DO  j = 1, nvex
-             a(ix) = a(ix+inc)
+             a(ix)     = a(ix+inc)
              a(ix+inc) = 0.0_wp
-             ix = ix + jump
+             ix        = ix + jump
           ENDDO
           IF ( MOD(n,2) /= 1 )  THEN
@@ -313,5 +331,5 @@
              DO  j = 1, nvex
                 a(iz) = 0.0_wp
-                iz = iz + jump
+                iz    = iz + jump
              ENDDO
           ENDIF
@@ -325,12 +343,11 @@
  END SUBROUTINE fft991cy
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Performs one pass through data as part of 
-!> multiple real fft (fourier analysis) routine.
-!> 
+!> Performs one pass through data as part of multiple real fft (fourier analysis) routine.
+!>
 !> Method:
-!> 
+!>
 !> a       is first real input vector
 !>         equivalence b(1) with a(ifac*la*inc1+1)
@@ -351,105 +368,105 @@
 !>         2 - ifac not catered for
 !>         3 - ifac only catered for if la=n/ifac
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE qpassm( a, b, c, d, trigs, inc1, inc2, inc3, inc4, lot, n, ifac, la, ierr )
 
     USE kinds
 
-    IMPLICIT NONE 
-
-    INTEGER(iwp) ::  ierr !< 
-    INTEGER(iwp) ::  ifac !< 
-    INTEGER(iwp) ::  inc1 !< 
-    INTEGER(iwp) ::  inc2 !< 
-    INTEGER(iwp) ::  inc3 !< 
-    INTEGER(iwp) ::  inc4 !< 
-    INTEGER(iwp) ::  la   !< 
-    INTEGER(iwp) ::  lot  !< 
-    INTEGER(iwp) ::  n    !< 
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  ierr  !<
+    INTEGER(iwp) ::  ifac  !<
+    INTEGER(iwp) ::  inc1  !<
+    INTEGER(iwp) ::  inc2  !<
+    INTEGER(iwp) ::  inc3  !<
+    INTEGER(iwp) ::  inc4  !<
+    INTEGER(iwp) ::  la    !<
+    INTEGER(iwp) ::  lot   !<
+    INTEGER(iwp) ::  n     !<
+
+    INTEGER(iwp) ::  i      !<
+    INTEGER(iwp) ::  ia     !<
+    INTEGER(iwp) ::  ib     !<
+    INTEGER(iwp) ::  ibase  !<
+    INTEGER(iwp) ::  ic     !<
+    INTEGER(iwp) ::  id     !<
+    INTEGER(iwp) ::  ie     !<
+    INTEGER(iwp) ::  if     !<
+    INTEGER(iwp) ::  ig     !<
+    INTEGER(iwp) ::  igo    !<
+    INTEGER(iwp) ::  ih     !<
+    INTEGER(iwp) ::  iink   !<
+    INTEGER(iwp) ::  ijk    !<
+    INTEGER(iwp) ::  ijump  !<
+    INTEGER(iwp) ::  j      !<
+    INTEGER(iwp) ::  ja     !<
+    INTEGER(iwp) ::  jb     !<
+    INTEGER(iwp) ::  jbase  !<
+    INTEGER(iwp) ::  jc     !<
+    INTEGER(iwp) ::  jd     !<
+    INTEGER(iwp) ::  je     !<
+    INTEGER(iwp) ::  jf     !<
+    INTEGER(iwp) ::  jink   !<
+    INTEGER(iwp) ::  k      !<
+    INTEGER(iwp) ::  kb     !<
+    INTEGER(iwp) ::  kc     !<
+    INTEGER(iwp) ::  kd     !<
+    INTEGER(iwp) ::  ke     !<
+    INTEGER(iwp) ::  kf     !<
+    INTEGER(iwp) ::  kstop  !<
+    INTEGER(iwp) ::  l      !<
+    INTEGER(iwp) ::  m      !<
 
 !
 !-- Arrays are dimensioned with n
-    REAL(wp) ::  a(*)     !< 
-    REAL(wp) ::  b(*)     !< 
-    REAL(wp) ::  c(*)     !< 
-    REAL(wp) ::  d(*)     !< 
-    REAL(wp) ::  trigs(*) !<
- 
-    REAL(wp) ::  a0     !< 
-    REAL(wp) ::  a1     !< 
-    REAL(wp) ::  a10    !< 
-    REAL(wp) ::  a11    !< 
-    REAL(wp) ::  a2     !< 
-    REAL(wp) ::  a20    !< 
-    REAL(wp) ::  a21    !< 
-    REAL(wp) ::  a3     !< 
-    REAL(wp) ::  a4     !< 
-    REAL(wp) ::  a5     !< 
-    REAL(wp) ::  a6     !< 
-    REAL(wp) ::  b0     !< 
-    REAL(wp) ::  b1     !< 
-    REAL(wp) ::  b10    !< 
-    REAL(wp) ::  b11    !< 
-    REAL(wp) ::  b2     !< 
-    REAL(wp) ::  b20    !< 
-    REAL(wp) ::  b21    !< 
-    REAL(wp) ::  b3     !< 
-    REAL(wp) ::  b4     !< 
-    REAL(wp) ::  b5     !< 
-    REAL(wp) ::  b6     !< 
-    REAL(wp) ::  c1     !< 
-    REAL(wp) ::  c2     !< 
-    REAL(wp) ::  c3     !< 
-    REAL(wp) ::  c4     !< 
-    REAL(wp) ::  c5     !< 
-    REAL(wp) ::  qrt5   !< 
-    REAL(wp) ::  s1     !< 
-    REAL(wp) ::  s2     !< 
-    REAL(wp) ::  s3     !< 
-    REAL(wp) ::  s4     !< 
-    REAL(wp) ::  s5     !< 
-    REAL(wp) ::  sin36  !< 
-    REAL(wp) ::  sin45  !< 
-    REAL(wp) ::  sin60  !< 
-    REAL(wp) ::  sin72  !< 
-    REAL(wp) ::  z      !< 
-    REAL(wp) ::  zqrt5  !< 
-    REAL(wp) ::  zsin36 !< 
-    REAL(wp) ::  zsin45 !< 
-    REAL(wp) ::  zsin60 !< 
-    REAL(wp) ::  zsin72 !< 
-
-    INTEGER(iwp) ::  i     !< 
-    INTEGER(iwp) ::  ia    !< 
-    INTEGER(iwp) ::  ib    !< 
-    INTEGER(iwp) ::  ibase !< 
-    INTEGER(iwp) ::  ic    !< 
-    INTEGER(iwp) ::  id    !< 
-    INTEGER(iwp) ::  ie    !< 
-    INTEGER(iwp) ::  if    !< 
-    INTEGER(iwp) ::  ig    !< 
-    INTEGER(iwp) ::  igo   !< 
-    INTEGER(iwp) ::  ih    !< 
-    INTEGER(iwp) ::  iink  !< 
-    INTEGER(iwp) ::  ijk   !< 
-    INTEGER(iwp) ::  ijump !< 
-    INTEGER(iwp) ::  j     !< 
-    INTEGER(iwp) ::  ja    !< 
-    INTEGER(iwp) ::  jb    !< 
-    INTEGER(iwp) ::  jbase !< 
-    INTEGER(iwp) ::  jc    !< 
-    INTEGER(iwp) ::  jd    !< 
-    INTEGER(iwp) ::  je    !< 
-    INTEGER(iwp) ::  jf    !< 
-    INTEGER(iwp) ::  jink  !< 
-    INTEGER(iwp) ::  k     !< 
-    INTEGER(iwp) ::  kb    !< 
-    INTEGER(iwp) ::  kc    !< 
-    INTEGER(iwp) ::  kd    !< 
-    INTEGER(iwp) ::  ke    !< 
-    INTEGER(iwp) ::  kf    !< 
-    INTEGER(iwp) ::  kstop !< 
-    INTEGER(iwp) ::  l     !< 
-    INTEGER(iwp) ::  m     !< 
+    REAL(wp) ::  a(*)      !<
+    REAL(wp) ::  b(*)      !<
+    REAL(wp) ::  c(*)      !<
+    REAL(wp) ::  d(*)      !<
+    REAL(wp) ::  trigs(*)  !<
+
+    REAL(wp) ::  a0      !<
+    REAL(wp) ::  a1      !<
+    REAL(wp) ::  a10     !<
+    REAL(wp) ::  a11     !<
+    REAL(wp) ::  a2      !<
+    REAL(wp) ::  a20     !<
+    REAL(wp) ::  a21     !<
+    REAL(wp) ::  a3      !<
+    REAL(wp) ::  a4      !<
+    REAL(wp) ::  a5      !<
+    REAL(wp) ::  a6      !<
+    REAL(wp) ::  b0      !<
+    REAL(wp) ::  b1      !<
+    REAL(wp) ::  b10     !<
+    REAL(wp) ::  b11     !<
+    REAL(wp) ::  b2      !<
+    REAL(wp) ::  b20     !<
+    REAL(wp) ::  b21     !<
+    REAL(wp) ::  b3      !<
+    REAL(wp) ::  b4      !<
+    REAL(wp) ::  b5      !<
+    REAL(wp) ::  b6      !<
+    REAL(wp) ::  c1      !<
+    REAL(wp) ::  c2      !<
+    REAL(wp) ::  c3      !<
+    REAL(wp) ::  c4      !<
+    REAL(wp) ::  c5      !<
+    REAL(wp) ::  qrt5    !<
+    REAL(wp) ::  s1      !<
+    REAL(wp) ::  s2      !<
+    REAL(wp) ::  s3      !<
+    REAL(wp) ::  s4      !<
+    REAL(wp) ::  s5      !<
+    REAL(wp) ::  sin36   !<
+    REAL(wp) ::  sin45   !<
+    REAL(wp) ::  sin60   !<
+    REAL(wp) ::  sin72   !<
+    REAL(wp) ::  z       !<
+    REAL(wp) ::  zqrt5   !<
+    REAL(wp) ::  zsin36  !<
+    REAL(wp) ::  zsin45  !<
+    REAL(wp) ::  zsin60  !<
+    REAL(wp) ::  zsin72  !<
 
 
@@ -468,6 +485,6 @@
     iink  = la * inc1
     jink  = la * inc2
-    ijump = (ifac-1) * iink
-    kstop = ( n-ifac ) / ( 2*ifac )
+    ijump = (ifac - 1) * iink
+    kstop = ( n - ifac ) / ( 2 * ifac )
 
     ibase = 0
@@ -490,5 +507,5 @@
           ib = ia + iink
           ja = 1
-          jb = ja + (2*m-la) * inc2
+          jb = ja + ( 2 * m - la ) * inc2
 
           IF ( la /= m )  THEN
@@ -510,7 +527,7 @@
              ENDDO
 
-             ja = ja + jink
-             jink = 2 * jink
-             jb = jb - jink
+             ja    = ja + jink
+             jink  = 2 * jink
+             jb    = jb - jink
              ibase = ibase + ijump
              ijump = 2 * ijump + iink
@@ -530,8 +547,8 @@
                       !OCL NOVREC
                       DO  ijk = 1, lot
-                         c(ja+j) = a(ia+i) + (c1*a(ib+i)+s1*b(ib+i))
-                         c(jb+j) = a(ia+i) - (c1*a(ib+i)+s1*b(ib+i))
-                         d(ja+j) = (c1*b(ib+i)-s1*a(ib+i)) + b(ia+i)
-                         d(jb+j) = (c1*b(ib+i)-s1*a(ib+i)) - b(ia+i)
+                         c(ja+j) = a(ia+i) + ( c1 * a(ib+i) + s1 * b(ib+i) )
+                         c(jb+j) = a(ia+i) - ( c1 * a(ib+i) + s1 * b(ib+i) )
+                         d(ja+j) = ( c1 * b(ib+i) - s1 * a(ib+i) ) + b(ia+i)
+                         d(jb+j) = ( c1 * b(ib+i) - s1 * a(ib+i) ) - b(ia+i)
                          i = i + inc3
                          j = j + inc4
@@ -562,5 +579,4 @@
                    i = i + inc3
                    j = j + inc4
-
                 ENDDO
                 ibase = ibase + inc1
@@ -570,5 +586,5 @@
           ELSE
 
-             z = 1.0_wp/REAL(n,KIND=wp)
+             z = 1.0_wp / REAL( n, KIND = wp )
              DO  l = 1, la
                 i = ibase
@@ -578,6 +594,6 @@
                 !OCL NOVREC
                 DO  ijk = 1, lot
-                   c(ja+j) = z*(a(ia+i)+a(ib+i))
-                   c(jb+j) = z*(a(ia+i)-a(ib+i))
+                   c(ja+j) = z * ( a(ia+i) + a(ib+i) )
+                   c(jb+j) = z * ( a(ia+i) - a(ib+i) )
                    i = i + inc3
                    j = j + inc4
@@ -597,5 +613,5 @@
           ic = ib + iink
           ja = 1
-          jb = ja + (2*m-la) * inc2
+          jb = ja + ( 2 * m - la) * inc2
           jc = jb
 
@@ -609,7 +625,7 @@
                 !OCL NOVREC
                 DO  ijk = 1, lot
-                   c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
-                   c(jb+j) = a(ia+i) - 0.5_wp*(a(ib+i)+a(ic+i))
-                   d(jb+j) = sin60*(a(ic+i)-a(ib+i))
+                   c(ja+j) = a(ia+i) + ( a(ib+i) + a(ic+i) )
+                   c(jb+j) = a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) )
+                   d(jb+j) = sin60 * ( a(ic+i) - a(ib+i) )
                    i = i + inc3
                    j = j + inc4
@@ -644,10 +660,12 @@
                       !OCL NOVREC
                       DO  ijk = 1, lot
-                         a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c2*a(ic+i)+s2*b(ic+i))
-                         b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c2*b(ic+i)-s2*a(ic+i))
-                         a2 = a(ia+i) - 0.5_wp*a1
-                         b2 = b(ia+i) - 0.5_wp*b1
-                         a3 = sin60*((c1*a(ib+i)+s1*b(ib+i))-(c2*a(ic+i)+s2*b(ic+i)))
-                         b3 = sin60*((c1*b(ib+i)-s1*a(ib+i))-(c2*b(ic+i)-s2*a(ic+i)))
+                         a1 = ( c1 * a(ib+i) + s1 * b(ib+i) ) + ( c2 * a(ic+i) + s2 * b(ic+i) )
+                         b1 = ( c1 * b(ib+i) - s1 * a(ib+i) ) + ( c2 * b(ic+i) - s2 * a(ic+i) )
+                         a2 = a(ia+i) - 0.5_wp * a1
+                         b2 = b(ia+i) - 0.5_wp * b1
+                         a3 = sin60 * ( ( c1 * a(ib+i) + s1 * b(ib+i) ) - ( c2 * a(ic+i) + s2      &
+                              * b(ic+i) ) )
+                         b3 = sin60 * ( ( c1 * b(ib+i) - s1 * a(ib+i) ) - ( c2 * b(ic+i) - s2      &
+                              * a(ic+i) ) )
                          c(ja+j) = a(ia+i) + a1
                          d(ja+j) = b(ia+i) + b1
@@ -655,5 +673,5 @@
                          d(jb+j) = b2 - a3
                          c(jc+j) = a2 - b3
-                         d(jc+j) = -(b2+a3)
+                         d(jc+j) = - ( b2 + a3 )
                          i = i + inc3
                          j = j + inc4
@@ -680,7 +698,7 @@
                 !OCL NOVREC
                 DO  ijk = 1, lot
-                   c(ja+j) = a(ia+i) + 0.5_wp*(a(ib+i)-a(ic+i))
-                   d(ja+j) = -sin60*(a(ib+i)+a(ic+i))
-                   c(jb+j) = a(ia+i) - (a(ib+i)-a(ic+i))
+                   c(ja+j) = a(ia+i) + 0.5_wp * ( a(ib+i) - a(ic+i) )
+                   d(ja+j) = -sin60 * ( a(ib+i) + a(ic+i) )
+                   c(jb+j) = a(ia+i) - ( a(ib+i) - a(ic+i) )
                    i = i + inc3
                    j = j + inc4
@@ -692,6 +710,6 @@
           ELSE
 
-             z = 1.0_wp / REAL( n, KIND=wp )
-             zsin60 = z*sin60
+             z = 1.0_wp / REAL( n, KIND = wp )
+             zsin60 = z * sin60
              DO  l = 1, la
                 i = ibase
@@ -701,7 +719,7 @@
                 !OCL NOVREC
                 DO  ijk = 1, lot
-                   c(ja+j) = z*(a(ia+i)+(a(ib+i)+a(ic+i)))
-                   c(jb+j) = z*(a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))
-                   d(jb+j) = zsin60*(a(ic+i)-a(ib+i))
+                   c(ja+j) = z * ( a(ia+i) + ( a (ib+i) + a(ic+i) ) )
+                   c(jb+j) = z * ( a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) ) )
+                   d(jb+j) = zsin60 * ( a(ic+i) -a(ib+i) )
                    i = i + inc3
                    j = j + inc4
@@ -722,6 +740,6 @@
           id = ic + iink
           ja = 1
-          jb = ja + (2*m-la) * inc2
-          jc = jb + 2*m*inc2
+          jb = ja + ( 2 * m - la ) * inc2
+          jc = jb + 2 * m * inc2
           jd = jb
 
@@ -736,6 +754,6 @@
                 !OCL NOVREC
                 DO  ijk = 1, lot
-                   c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))
-                   c(jc+j) = (a(ia+i)+a(ic+i)) - (a(ib+i)+a(id+i))
+                   c(ja+j) = ( a(ia+i) + a(ic+i) ) + ( a(ib+i) + a(id+i) )
+                   c(jc+j) = ( a(ia+i) + a(ic+i) ) - ( a(ib+i) + a(id+i) )
                    c(jb+j) = a(ia+i) - a(ic+i)
                    d(jb+j) = a(id+i) - a(ib+i)
@@ -747,9 +765,9 @@
              ENDDO
 
-             ja = ja + jink
-             jink = 2 * jink
-             jb = jb + jink
-             jc = jc - jink
-             jd = jd - jink
+             ja    = ja + jink
+             jink  = 2 * jink
+             jb    = jb + jink
+             jc    = jc - jink
+             jd    = jd - jink
              ibase = ibase + ijump
              ijump = 2 * ijump + iink
@@ -775,12 +793,12 @@
                       !OCL NOVREC
                       DO  ijk = 1, lot
-                         a0 = a(ia+i) + (c2*a(ic+i)+s2*b(ic+i))
-                         a2 = a(ia+i) - (c2*a(ic+i)+s2*b(ic+i))
-                         a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c3*a(id+i)+s3*b(id+i))
-                         a3 = (c1*a(ib+i)+s1*b(ib+i)) - (c3*a(id+i)+s3*b(id+i))
-                         b0 = b(ia+i) + (c2*b(ic+i)-s2*a(ic+i))
-                         b2 = b(ia+i) - (c2*b(ic+i)-s2*a(ic+i))
-                         b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c3*b(id+i)-s3*a(id+i))
-                         b3 = (c1*b(ib+i)-s1*a(ib+i)) - (c3*b(id+i)-s3*a(id+i))
+                         a0 = a(ia+i) + ( c2 * a(ic+i) + s2 * b(ic+i) )
+                         a2 = a(ia+i) - ( c2 * a(ic+i) + s2 * b(ic+i) )
+                         a1 = ( c1 * a(ib+i) + s1 * b(ib+i) ) + ( c3 * a(id+i) + s3 * b(id+i) )
+                         a3 = ( c1 * a(ib+i) + s1 * b(ib+i) ) - ( c3 * a(id+i) + s3 * b(id+i) )
+                         b0 = b(ia+i) + ( c2 * b(ic+i) - s2 * a(ic+i) )
+                         b2 = b(ia+i) - ( c2 * b(ic+i) - s2 * a(ic+i) )
+                         b1 = ( c1 * b(ib+i) - s1 * a(ib+i) ) + ( c3 * b(id+i) - s3 * a(id+i) )
+                         b3 = ( c1 * b(ib+i) - s1 * a(ib+i) ) - ( c3 * b(id+i) - s3 * a(id+i) )
                          c(ja+j) = a0 + a1
                          c(jc+j) = a0 - a1
@@ -790,5 +808,5 @@
                          c(jd+j) = a2 - b3
                          d(jb+j) = b2 - a3
-                         d(jd+j) = -(b2+a3)
+                         d(jd+j) = - ( b2 + a3 )
                          i = i + inc3
                          j = j + inc4
@@ -817,8 +835,8 @@
                 !OCL NOVREC
                 DO  ijk = 1, lot
-                   c(ja+j) = a(ia+i) + sin45*(a(ib+i)-a(id+i))
-                   c(jb+j) = a(ia+i) - sin45*(a(ib+i)-a(id+i))
-                   d(ja+j) = -a(ic+i) - sin45*(a(ib+i)+a(id+i))
-                   d(jb+j) = a(ic+i) - sin45*(a(ib+i)+a(id+i))
+                   c(ja+j) = a(ia+i) + sin45 * ( a(ib+i) - a(id+i) )
+                   c(jb+j) = a(ia+i) - sin45 * ( a(ib+i) - a(id+i) )
+                   d(ja+j) = - a(ic+i) - sin45 * ( a(ib+i) + a(id+i) )
+                   d(jb+j) = a(ic+i) - sin45 * ( a(ib+i) + a(id+i) )
                    i = i + inc3
                    j = j + inc4
@@ -830,5 +848,5 @@
           ELSE
 
-             z = 1.0_wp / REAL( n, KIND=wp )
+             z = 1.0_wp / REAL( n, KIND = wp )
              DO  l = 1, la
                 i = ibase
@@ -838,8 +856,8 @@
                 !OCL NOVREC
                 DO  ijk = 1, lot
-                   c(ja+j) = z*((a(ia+i)+a(ic+i))+(a(ib+i)+a(id+i)))
-                   c(jc+j) = z*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i)))
-                   c(jb+j) = z*(a(ia+i)-a(ic+i))
-                   d(jb+j) = z*(a(id+i)-a(ib+i))
+                   c(ja+j) = z * ( ( a(ia+i) + a(ic+i) ) + ( a(ib+i) + a(id+i) ) )
+                   c(jc+j) = z * ( ( a(ia+i) + a(ic+i) ) - (a (ib+i) + a(id+i) ) )
+                   c(jb+j) = z * ( a(ia+i) - a(ic+i) )
+                   d(jb+j) = z * ( a(id+i) - a(ib+i) )
                    i = i + inc3
                    j = j + inc4
@@ -861,6 +879,6 @@
           ie = id + iink
           ja = 1
-          jb = ja + (2*m-la) * inc2
-          jc = jb + 2*m*inc2
+          jb = ja + ( 2 * m - la ) * inc2
+          jc = jb + 2 * m * inc2
           jd = jc
           je = jb
@@ -879,11 +897,11 @@
                    a2 = a(ic+i) + a(id+i)
                    a4 = a(ic+i) - a(id+i)
-                   a5 = a(ia+i) - 0.25_wp*(a1+a2)
-                   a6 = qrt5*(a1-a2)
-                   c(ja+j) = a(ia+i) + (a1+a2)
+                   a5 = a(ia+i) - 0.25_wp * ( a1 + a2 )
+                   a6 = qrt5 * ( a1 - a2 )
+                   c(ja+j) = a(ia+i) + ( a1 + a2 )
                    c(jb+j) = a5 + a6
                    c(jc+j) = a5 - a6
-                   d(jb+j) = -sin72*a3 - sin36*a4
-                   d(jc+j) = -sin36*a3 + sin72*a4
+                   d(jb+j) = - sin72 * a3 - sin36 * a4
+                   d(jc+j) = - sin36 * a3 + sin72 * a4
                    i = i + inc3
                    j = j + inc4
@@ -925,34 +943,34 @@
                       !OCL NOVREC
                       DO  ijk = 1, lot
-                         a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c4*a(ie+i)+s4*b(ie+i))
-                         a3 = (c1*a(ib+i)+s1*b(ib+i)) - (c4*a(ie+i)+s4*b(ie+i))
-                         a2 = (c2*a(ic+i)+s2*b(ic+i)) + (c3*a(id+i)+s3*b(id+i))
-                         a4 = (c2*a(ic+i)+s2*b(ic+i)) - (c3*a(id+i)+s3*b(id+i))
-                         b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c4*b(ie+i)-s4*a(ie+i))
-                         b3 = (c1*b(ib+i)-s1*a(ib+i)) - (c4*b(ie+i)-s4*a(ie+i))
-                         b2 = (c2*b(ic+i)-s2*a(ic+i)) + (c3*b(id+i)-s3*a(id+i))
-                         b4 = (c2*b(ic+i)-s2*a(ic+i)) - (c3*b(id+i)-s3*a(id+i))
-                         a5 = a(ia+i) - 0.25_wp*(a1+a2)
-                         a6 = qrt5*(a1-a2)
-                         b5 = b(ia+i) - 0.25_wp*(b1+b2)
-                         b6 = qrt5*(b1-b2)
+                         a1  = (c1 * a(ib+i) + s1 * b(ib+i) ) + ( c4 * a(ie+i) + s4 * b(ie+i) )
+                         a3  = (c1 * a(ib+i) + s1 * b(ib+i) ) - ( c4 * a(ie+i) + s4 * b(ie+i) )
+                         a2  = (c2 * a(ic+i) + s2 * b(ic+i) ) + ( c3 * a(id+i) + s3 * b(id+i) )
+                         a4  = (c2 * a(ic+i) + s2 * b(ic+i) ) - ( c3 * a(id+i) + s3 * b(id+i) )
+                         b1  = (c1 * b(ib+i) - s1 * a(ib+i) ) + ( c4 * b(ie+i) - s4 * a(ie+i) )
+                         b3  = (c1 * b(ib+i) - s1 * a(ib+i) ) - ( c4 * b(ie+i) - s4 * a(ie+i) )
+                         b2  = (c2 * b(ic+i) - s2 * a(ic+i) ) + ( c3 * b(id+i) - s3 * a(id+i) )
+                         b4  = (c2 * b(ic+i) - s2 * a(ic+i) ) - ( c3 * b(id+i) - s3 * a(id+i) )
+                         a5  = a(ia+i) - 0.25_wp * ( a1 + a2 )
+                         a6  = qrt5 * ( a1 - a2 )
+                         b5  = b(ia+i) - 0.25_wp * ( b1 + b2 )
+                         b6  = qrt5 * ( b1 - b2 )
                          a10 = a5 + a6
                          a20 = a5 - a6
                          b10 = b5 + b6
                          b20 = b5 - b6
-                         a11 = sin72*b3 + sin36*b4
-                         a21 = sin36*b3 - sin72*b4
-                         b11 = sin72*a3 + sin36*a4
-                         b21 = sin36*a3 - sin72*a4
-                         c(ja+j) = a(ia+i) + (a1+a2)
+                         a11 = sin72 * b3 + sin36 * b4
+                         a21 = sin36 * b3 - sin72 * b4
+                         b11 = sin72 * a3 + sin36 * a4
+                         b21 = sin36 * a3 - sin72 * a4
+                         c(ja+j) = a(ia+i) + ( a1 + a2 )
                          c(jb+j) = a10 + a11
                          c(je+j) = a10 - a11
                          c(jc+j) = a20 + a21
                          c(jd+j) = a20 - a21
-                         d(ja+j) = b(ia+i) + (b1+b2)
+                         d(ja+j) = b(ia+i) + ( b1 + b2 )
                          d(jb+j) = b10 - b11
-                         d(je+j) = -(b10+b11)
+                         d(je+j) = - ( b10 + b11 )
                          d(jc+j) = b20 - b21
-                         d(jd+j) = -(b20+b21)
+                         d(jd+j) = - ( b20 + b21 )
                          i = i + inc3
                          j = j + inc4
@@ -985,11 +1003,11 @@
                    a2 = a(ic+i) + a(id+i)
                    a4 = a(ic+i) - a(id+i)
-                   a5 = a(ia+i) + 0.25_wp*(a3-a4)
-                   a6 = qrt5*(a3+a4)
+                   a5 = a(ia+i) + 0.25_wp * ( a3 - a4 )
+                   a6 = qrt5 * ( a3 + a4 )
                    c(ja+j) = a5 + a6
                    c(jb+j) = a5 - a6
-                   c(jc+j) = a(ia+i) - (a3-a4)
-                   d(ja+j) = -sin36*a1 - sin72*a2
-                   d(jb+j) = -sin72*a1 + sin36*a2
+                   c(jc+j) = a(ia+i) - ( a3 - a4 )
+                   d(ja+j) = - sin36 * a1 - sin72 * a2
+                   d(jb+j) = - sin72 * a1 + sin36 * a2
                    i = i + inc3
                    j = j + inc4
@@ -1001,5 +1019,5 @@
           ELSE
 
-             z = 1.0_wp / REAL( n, KIND=wp )
+             z = 1.0_wp / REAL( n, KIND = wp )
              zqrt5  = z * qrt5
              zsin36 = z * sin36
@@ -1016,11 +1034,11 @@
                    a2 = a(ic+i) + a(id+i)
                    a4 = a(ic+i) - a(id+i)
-                   a5 = z*(a(ia+i)-0.25_wp*(a1+a2))
-                   a6 = zqrt5*(a1-a2)
-                   c(ja+j) = z*(a(ia+i)+(a1+a2))
+                   a5 = z * ( a (ia+i) - 0.25_wp * ( a1 + a2 ) )
+                   a6 = zqrt5 * ( a1 - a2 )
+                   c(ja+j) = z * ( a(ia+i) + (a1+a2) )
                    c(jb+j) = a5 + a6
                    c(jc+j) = a5 - a6
-                   d(jb+j) = -zsin72*a3 - zsin36*a4
-                   d(jc+j) = -zsin36*a3 + zsin72*a4
+                   d(jb+j) = - zsin72 * a3 - zsin36 * a4
+                   d(jc+j) = - zsin36 * a3 + zsin72 * a4
                    i = i + inc3
                    j = j + inc4
@@ -1043,7 +1061,7 @@
           if = ie + iink
           ja = 1
-          jb = ja + (2*m-la) * inc2
-          jc = jb + 2*m*inc2
-          jd = jc + 2*m*inc2
+          jb = ja + ( 2 * m - la ) * inc2
+          jc = jb + 2 * m * inc2
+          jd = jc + 2 * m * inc2
           je = jc
           jf = jb
@@ -1058,12 +1076,12 @@
                 !OCL NOVREC
                 DO  ijk = 1, lot
-                   a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
-                   c(ja+j) = (a(ia+i)+a(id+i)) + a11
-                   c(jc+j) = (a(ia+i)+a(id+i)-0.5_wp*a11)
-                   d(jc+j) = sin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
-                   a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
-                   c(jb+j) = (a(ia+i)-a(id+i)) - 0.5_wp*a11
-                   d(jb+j) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
-                   c(jd+j) = (a(ia+i)-a(id+i)) + a11
+                   a11 = ( a (ic+i) + a(if+i) ) + ( a(ib+i) + a(ie+i) )
+                   c(ja+j) = ( a(ia+i) + a(id+i) ) + a11
+                   c(jc+j) = ( a(ia+i) + a(id+i) - 0.5_wp * a11 )
+                   d(jc+j) = sin60 * ( ( a(ic+i) + a(if+i) ) - ( a(ib+i) + a(ie+i) ) )
+                   a11 = ( a (ic+i) - a(if+i) ) + ( a(ie+i) - a(ib+i) )
+                   c(jb+j) = ( a(ia+i) - a(id+i) ) - 0.5_wp * a11
+                   d(jb+j) = sin60 * ( ( a(ie+i) - a(ib+i) ) - ( a(ic+i) - a(if+i) ) )
+                   c(jd+j) = ( a (ia+i) - a(id+i) ) + a11
                    i = i + inc3
                    j = j + inc4
@@ -1072,11 +1090,11 @@
                 jbase = jbase + inc2
              ENDDO
-             ja = ja + jink
-             jink = 2 * jink
-             jb = jb + jink
-             jc = jc + jink
-             jd = jd - jink
-             je = je - jink
-             jf = jf - jink
+             ja    = ja + jink
+             jink  = 2 * jink
+             jb    = jb + jink
+             jc    = jc + jink
+             jd    = jd - jink
+             je    = je - jink
+             jf    = jf - jink
              ibase = ibase + ijump
              ijump = 2 * ijump + iink
@@ -1108,36 +1126,36 @@
                       !OCL NOVREC
                       DO  ijk = 1, lot
-                         a1 = c1*a(ib+i) + s1*b(ib+i)
-                         b1 = c1*b(ib+i) - s1*a(ib+i)
-                         a2 = c2*a(ic+i) + s2*b(ic+i)
-                         b2 = c2*b(ic+i) - s2*a(ic+i)
-                         a3 = c3*a(id+i) + s3*b(id+i)
-                         b3 = c3*b(id+i) - s3*a(id+i)
-                         a4 = c4*a(ie+i) + s4*b(ie+i)
-                         b4 = c4*b(ie+i) - s4*a(ie+i)
-                         a5 = c5*a(if+i) + s5*b(if+i)
-                         b5 = c5*b(if+i) - s5*a(if+i)
-                         a11 = (a2+a5) + (a1+a4)
-                         a20 = (a(ia+i)+a3) - 0.5_wp*a11
-                         a21 = sin60*((a2+a5)-(a1+a4))
-                         b11 = (b2+b5) + (b1+b4)
-                         b20 = (b(ia+i)+b3) - 0.5_wp*b11
-                         b21 = sin60*((b2+b5)-(b1+b4))
-                         c(ja+j) = (a(ia+i)+a3) + a11
-                         d(ja+j) = (b(ia+i)+b3) + b11
+                         a1  = c1 * a(ib+i) + s1 * b(ib+i)
+                         b1  = c1 * b(ib+i) - s1 * a(ib+i)
+                         a2  = c2 * a(ic+i) + s2 * b(ic+i)
+                         b2  = c2 * b(ic+i) - s2 * a(ic+i)
+                         a3  = c3 * a(id+i) + s3 * b(id+i)
+                         b3  = c3 * b(id+i) - s3 * a(id+i)
+                         a4  = c4 * a(ie+i) + s4 * b(ie+i)
+                         b4  = c4 * b(ie+i) - s4 * a(ie+i)
+                         a5  = c5 * a(if+i) + s5 * b(if+i)
+                         b5  = c5 * b(if+i) - s5 * a(if+i)
+                         a11 = ( a2 + a5 ) + ( a1 + a4 )
+                         a20 = ( a(ia+i) + a3 ) - 0.5_wp * a11
+                         a21 = sin60 * ( ( a2 + a5 ) - ( a1 + a4 ) )
+                         b11 = ( b2 + b5 ) + ( b1 + b4 )
+                         b20 = ( b(ia+i) + b3 ) - 0.5_wp * b11
+                         b21 = sin60 * ( ( b2 + b5 ) - ( b1 + b4 ) )
+                         c(ja+j) = ( a(ia+i) + a3 ) + a11
+                         d(ja+j) = ( b(ia+i) + b3 ) + b11
                          c(jc+j) = a20 - b21
                          d(jc+j) = a21 + b20
                          c(je+j) = a20 + b21
                          d(je+j) = a21 - b20
-                         a11 = (a2-a5) + (a4-a1)
-                         a20 = (a(ia+i)-a3) - 0.5_wp*a11
-                         a21 = sin60*((a4-a1)-(a2-a5))
-                         b11 = (b5-b2) - (b4-b1)
-                         b20 = (b3-b(ia+i)) - 0.5_wp*b11
-                         b21 = sin60*((b5-b2)+(b4-b1))
+                         a11 = ( a2 - a5 ) + ( a4 - a1 )
+                         a20 = ( a(ia+i) - a3 ) - 0.5_wp * a11
+                         a21 = sin60 * ( ( a4 - a1 ) - ( a2 - a5 ) )
+                         b11 = ( b5 - b2 ) - ( b4 - b1 )
+                         b20 = ( b3 - b(ia+i) ) - 0.5_wp * b11
+                         b21 = sin60 * ( ( b5 - b2 ) + ( b4 - b1 ) )
                          c(jb+j) = a20 - b21
                          d(jb+j) = a21 - b20
-                         c(jd+j) = a11 + (a(ia+i)-a3)
-                         d(jd+j) = b11 + (b3-b(ia+i))
+                         c(jd+j) = a11 + ( a(ia+i) - a3 )
+                         d(jd+j) = b11 + ( b3 - b(ia+i) )
                          c(jf+j) = a20 + b21
                          d(jf+j) = a21 + b20
@@ -1169,10 +1187,14 @@
                 !OCL NOVREC
                 DO  ijk = 1, lot
-                   c(ja+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) + sin60*(a(ib+i)-a(if+i))
-                   d(ja+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) - sin60*(a(ic+i)+a(ie+i))
-                   c(jb+j) = a(ia+i) - (a(ic+i)-a(ie+i))
-                   d(jb+j) = a(id+i) - (a(ib+i)+a(if+i))
-                   c(jc+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) - sin60*(a(ib+i)-a(if+i))
-                   d(jc+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) + sin60*(a(ic+i)+a(ie+i))
+                   c(ja+j) = ( a (ia+i) + 0.5_wp * ( a(ic+i) - a(ie+i) ) ) + sin60 *               &
+                             ( a(ib+i) - a(if+i) )
+                   d(ja+j) = - ( a(id+i) + 0.5_wp * ( a(ib+i) + a(if+i) ) ) - sin60 * ( a(ic+i) +  &
+                             a(ie+i) )
+                   c(jb+j) = a(ia+i) - ( a(ic+i) - a(ie+i) )
+                   d(jb+j) = a(id+i) - ( a(ib+i) + a(if+i) )
+                   c(jc+j) = ( a(ia+i) + 0.5_wp * ( a(ic+i) - a(ie+i) ) ) - sin60 * ( a(ib+i) -    &
+                             a(if+i) )
+                   d(jc+j) = - ( a(id+i) + 0.5_wp * ( a(ib+i) + a(if+i) ) ) + sin60 * ( a(ic+i) +  &
+                             a(ie+i) )
                    i = i + inc3
                    j = j + inc4
@@ -1184,6 +1206,6 @@
           ELSE
 
-             z = 1.0_wp/REAL(n,KIND=wp)
-             zsin60 = z*sin60
+             z = 1.0_wp / REAL( n, KIND = wp )
+             zsin60 = z * sin60
              DO  l = 1, la
                 i = ibase
@@ -1193,12 +1215,12 @@
                 !OCL NOVREC
                 DO  ijk = 1, lot
-                   a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
-                   c(ja+j) = z*((a(ia+i)+a(id+i))+a11)
-                   c(jc+j) = z*((a(ia+i)+a(id+i))-0.5_wp*a11)
-                   d(jc+j) = zsin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
-                   a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
-                   c(jb+j) = z*((a(ia+i)-a(id+i))-0.5_wp*a11)
-                   d(jb+j) = zsin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
-                   c(jd+j) = z*((a(ia+i)-a(id+i))+a11)
+                   a11 = ( a(ic+i) + a(if+i) ) + ( a(ib+i) + a(ie+i) )
+                   c(ja+j) = z * ( ( a(ia+i) + a(id+i) ) + a11 )
+                   c(jc+j) = z * ( ( a(ia+i) + a(id+i) ) - 0.5_wp * a11 )
+                   d(jc+j) = zsin60 * ( ( a(ic+i) + a(if+i) ) - ( a(ib+i) + a(ie+i) ) )
+                   a11 = ( a(ic+i) - a(if+i) ) + ( a(ie+i) - a(ib+i) )
+                   c(jb+j) = z * ( ( a(ia+i) - a(id+i) ) - 0.5_wp * a11 )
+                   d(jb+j) = zsin60 * ( ( a(ie+i) - a(ib+i) ) - ( a(ic+i) - a(if+i) ) )
+                   c(jd+j) = z * ( ( a(ia+i) - a(id+i) ) + a11 )
                    i = i + inc3
                    j = j + inc4
@@ -1229,8 +1251,8 @@
           ja = 1
           jb = ja + la * inc2
-          jc = jb + 2*m*inc2
-          jd = jc + 2*m*inc2
-          je = jd + 2*m*inc2
-          z = 1.0_wp / REAL( n, KIND=wp )
+          jc = jb + 2 * m * inc2
+          jd = jc + 2 * m * inc2
+          je = jd + 2 * m * inc2
+          z = 1.0_wp / REAL( n, KIND = wp )
           zsin45 = z * SQRT( 0.5_wp )
 
@@ -1243,12 +1265,18 @@
              !OCL NOVREC
              DO  ijk = 1, lot
-                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))))
-                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))))
-                c(jc+j) = z*((a(ia+i)+a(ie+i))-(a(ic+i)+a(ig+i)))
-                d(jc+j) = z*((a(id+i)+a(ih+i))-(a(ib+i)+a(if+i)))
-                c(jb+j) = z*(a(ia+i)-a(ie+i)) + zsin45*((a(ih+i)-a(id+i))-(a(if+i)-a(ib+i)))
-                c(jd+j) = z*(a(ia+i)-a(ie+i)) - zsin45*((a(ih+i)-a(id+i))-(a(if+i)-a(ib+i)))
-                d(jb+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) + z*(a(ig+i)-a(ic+i))
-                d(jd+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) - z*(a(ig+i)-a(ic+i))
+                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) ) ) )
+                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) ) ) )
+                c(jc+j) = z * ( ( a(ia+i) + a(ie+i) ) - ( a(ic+i) + a(ig+i) ) )
+                d(jc+j) = z * ( ( a(id+i) + a(ih+i) ) - ( a(ib+i) + a(if+i) ) )
+                c(jb+j) = z * ( a(ia+i) - a(ie+i) ) + zsin45 * ( ( a(ih+i) - a(id+i) ) - ( a(if+i) &
+                          - a(ib+i) ) )
+                c(jd+j) = z * ( a(ia+i) - a(ie+i) ) - zsin45 * ( ( a(ih+i) - a(id+i) ) - ( a(if+i) &
+                          - a(ib+i) ) )
+                d(jb+j) = zsin45 * ( ( a(ih+i) - a(id+i) ) + ( a(if+i) - a(ib+i) ) ) + z *         &
+                          ( a(ig+i) - a(ic+i) )
+                d(jd+j) = zsin45 * ( ( a(ih+i) - a(id+i) ) + ( a(if+i) - a(ib+i) ) ) - z *         &
+                          ( a(ig+i) - a(ic+i) )
                 i = i + inc3
                 j = j + inc4
@@ -1262,9 +1290,9 @@
  END SUBROUTINE qpassm
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Same as qpassm, but for backward fft
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE rpassm( a, b, c, d, trigs, inc1, inc2, inc3, inc4, lot, n, ifac, la, ierr )
 
@@ -1273,83 +1301,83 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  ierr !< 
-    INTEGER(iwp) ::  ifac !< 
-    INTEGER(iwp) ::  inc1 !< 
-    INTEGER(iwp) ::  inc2 !< 
-    INTEGER(iwp) ::  inc3 !< 
-    INTEGER(iwp) ::  inc4 !< 
-    INTEGER(iwp) ::  la   !< 
-    INTEGER(iwp) ::  lot  !< 
-    INTEGER(iwp) ::  n    !< 
+    INTEGER(iwp) ::  ierr  !<
+    INTEGER(iwp) ::  ifac  !<
+    INTEGER(iwp) ::  inc1  !<
+    INTEGER(iwp) ::  inc2  !<
+    INTEGER(iwp) ::  inc3  !<
+    INTEGER(iwp) ::  inc4  !<
+    INTEGER(iwp) ::  la    !<
+    INTEGER(iwp) ::  lot   !<
+    INTEGER(iwp) ::  n     !<
 !
 !-- Arrays are dimensioned with n
-    REAL(wp) ::  a(*)     !< 
-    REAL(wp) ::  b(*)     !< 
-    REAL(wp) ::  c(*)     !< 
-    REAL(wp) ::  d(*)     !< 
-    REAL(wp) ::  trigs(*) !< 
-
-    REAL(wp) ::  c1     !< 
-    REAL(wp) ::  c2     !< 
-    REAL(wp) ::  c3     !< 
-    REAL(wp) ::  c4     !< 
-    REAL(wp) ::  c5     !< 
-    REAL(wp) ::  qqrt5  !< 
-    REAL(wp) ::  qrt5   !< 
-    REAL(wp) ::  s1     !< 
-    REAL(wp) ::  s2     !< 
-    REAL(wp) ::  s3     !< 
-    REAL(wp) ::  s4     !< 
-    REAL(wp) ::  s5     !< 
-    REAL(wp) ::  sin36  !< 
-    REAL(wp) ::  sin45  !< 
-    REAL(wp) ::  sin60  !< 
-    REAL(wp) ::  sin72  !< 
-    REAL(wp) ::  ssin36 !< 
-    REAL(wp) ::  ssin45 !< 
-    REAL(wp) ::  ssin60 !< 
-    REAL(wp) ::  ssin72 !< 
-
-    INTEGER(iwp) ::  i     !< 
-    INTEGER(iwp) ::  ia    !< 
-    INTEGER(iwp) ::  ib    !< 
-    INTEGER(iwp) ::  ibase !< 
-    INTEGER(iwp) ::  ic    !< 
-    INTEGER(iwp) ::  id    !< 
-    INTEGER(iwp) ::  ie    !< 
-    INTEGER(iwp) ::  if    !< 
-    INTEGER(iwp) ::  igo   !< 
-    INTEGER(iwp) ::  iink  !< 
-    INTEGER(iwp) ::  ijk   !< 
-    INTEGER(iwp) ::  j     !< 
-    INTEGER(iwp) ::  ja    !< 
-    INTEGER(iwp) ::  jb    !< 
-    INTEGER(iwp) ::  jbase !< 
-    INTEGER(iwp) ::  jc    !< 
-    INTEGER(iwp) ::  jd    !< 
-    INTEGER(iwp) ::  je    !< 
-    INTEGER(iwp) ::  jf    !< 
-    INTEGER(iwp) ::  jg    !< 
-    INTEGER(iwp) ::  jh    !< 
-    INTEGER(iwp) ::  jink  !< 
-    INTEGER(iwp) ::  jump  !< 
-    INTEGER(iwp) ::  k     !< 
-    INTEGER(iwp) ::  kb    !< 
-    INTEGER(iwp) ::  kc    !< 
-    INTEGER(iwp) ::  kd    !< 
-    INTEGER(iwp) ::  ke    !< 
-    INTEGER(iwp) ::  kf    !< 
-    INTEGER(iwp) ::  kstop !< 
-    INTEGER(iwp) ::  l     !< 
-    INTEGER(iwp) ::  m     !< 
-
-    REAL(wp) ::  a10(nfft) !< 
-    REAL(wp) ::  a11(nfft) !< 
-    REAL(wp) ::  a20(nfft) !< 
-    REAL(wp) ::  a21(nfft) !< 
-    REAL(wp) ::  b10(nfft) !< 
-    REAL(wp) ::  b11(nfft) !< 
-    REAL(wp) ::  b20(nfft) !< 
-    REAL(wp) ::  b21(nfft) !< 
+    INTEGER(iwp) ::  i      !<
+    INTEGER(iwp) ::  ia     !<
+    INTEGER(iwp) ::  ib     !<
+    INTEGER(iwp) ::  ibase  !<
+    INTEGER(iwp) ::  ic     !<
+    INTEGER(iwp) ::  id     !<
+    INTEGER(iwp) ::  ie     !<
+    INTEGER(iwp) ::  if     !<
+    INTEGER(iwp) ::  igo    !<
+    INTEGER(iwp) ::  iink   !<
+    INTEGER(iwp) ::  ijk    !<
+    INTEGER(iwp) ::  j      !<
+    INTEGER(iwp) ::  ja     !<
+    INTEGER(iwp) ::  jb     !<
+    INTEGER(iwp) ::  jbase  !<
+    INTEGER(iwp) ::  jc     !<
+    INTEGER(iwp) ::  jd     !<
+    INTEGER(iwp) ::  je     !<
+    INTEGER(iwp) ::  jf     !<
+    INTEGER(iwp) ::  jg     !<
+    INTEGER(iwp) ::  jh     !<
+    INTEGER(iwp) ::  jink   !<
+    INTEGER(iwp) ::  jump   !<
+    INTEGER(iwp) ::  k      !<
+    INTEGER(iwp) ::  kb     !<
+    INTEGER(iwp) ::  kc     !<
+    INTEGER(iwp) ::  kd     !<
+    INTEGER(iwp) ::  ke     !<
+    INTEGER(iwp) ::  kf     !<
+    INTEGER(iwp) ::  kstop  !<
+    INTEGER(iwp) ::  l      !<
+    INTEGER(iwp) ::  m      !<
+
+    REAL(wp) ::  a(*)      !<
+    REAL(wp) ::  b(*)      !<
+    REAL(wp) ::  c(*)      !<
+    REAL(wp) ::  d(*)      !<
+    REAL(wp) ::  trigs(*)  !<
+
+    REAL(wp) ::  c1      !<
+    REAL(wp) ::  c2      !<
+    REAL(wp) ::  c3      !<
+    REAL(wp) ::  c4      !<
+    REAL(wp) ::  c5      !<
+    REAL(wp) ::  qqrt5   !<
+    REAL(wp) ::  qrt5    !<
+    REAL(wp) ::  s1      !<
+    REAL(wp) ::  s2      !<
+    REAL(wp) ::  s3      !<
+    REAL(wp) ::  s4      !<
+    REAL(wp) ::  s5      !<
+    REAL(wp) ::  sin36   !<
+    REAL(wp) ::  sin45   !<
+    REAL(wp) ::  sin60   !<
+    REAL(wp) ::  sin72   !<
+    REAL(wp) ::  ssin36  !<
+    REAL(wp) ::  ssin45  !<
+    REAL(wp) ::  ssin60  !<
+    REAL(wp) ::  ssin72  !<
+
+    REAL(wp) ::  a10(nfft)  !<
+    REAL(wp) ::  a11(nfft)  !<
+    REAL(wp) ::  a20(nfft)  !<
+    REAL(wp) ::  a21(nfft)  !<
+    REAL(wp) ::  b10(nfft)  !<
+    REAL(wp) ::  b11(nfft)  !<
+    REAL(wp) ::  b20(nfft)  !<
+    REAL(wp) ::  b21(nfft)  !<
 
 
@@ -1363,6 +1391,6 @@
     iink = la * inc1
     jink = la * inc2
-    jump = (ifac-1) * jink
-    kstop = (n-ifac) / (2*ifac)
+    jump = ( ifac - 1 ) * jink
+    kstop = ( n - ifac ) / ( 2 * ifac )
 
     IF ( lot > nfft )  THEN
@@ -1386,5 +1414,5 @@
 
           ia = 1
-          ib = ia + (2*m-la) * inc1
+          ib = ia + ( 2 * m - la ) * inc1
           ja = 1
           jb = ja + jink
@@ -1409,5 +1437,5 @@
 
              ia = ia + iink
-             iink = 2 * iink
+             iink  = 2 * iink
              ib = ib - iink
              ibase = 0
@@ -1431,6 +1459,6 @@
                          c(ja+j) = a(ia+i) + a(ib+i)
                          d(ja+j) = b(ia+i) - b(ib+i)
-                         c(jb+j) = c1*(a(ia+i)-a(ib+i)) - s1*(b(ia+i)+b(ib+i))
-                         d(jb+j) = s1*(a(ia+i)-a(ib+i)) + c1*(b(ia+i)+b(ib+i))
+                         c(jb+j) = c1 * ( a(ia+i) - a(ib+i) ) - s1 * ( b(ia+i) + b(ib+i) )
+                         d(jb+j) = s1 * ( a(ia+i) - a(ib+i) ) + c1 * ( b(ia+i) + b(ib+i) )
                          i = i + inc3
                          j = j + inc4
@@ -1457,5 +1485,5 @@
                 DO  ijk = 1, lot
                    c(ja+j) = a(ia+i)
-                   c(jb+j) = -b(ia+i)
+                   c(jb+j) = - b(ia+i)
                    i = i + inc3
                    j = j + inc4
@@ -1474,6 +1502,6 @@
                 !OCL NOVREC
                 DO  ijk = 1, lot
-                   c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i))
-                   c(jb+j) = 2.0_wp*(a(ia+i)-a(ib+i))
+                   c(ja+j) = 2.0_wp * ( a(ia+i) + a(ib+i) )
+                   c(jb+j) = 2.0_wp * ( a(ia+i) - a(ib+i) )
                    i = i + inc3
                    j = j + inc4
@@ -1490,5 +1518,5 @@
 
           ia = 1
-          ib = ia + (2*m-la) * inc1
+          ib = ia + ( 2 * m - la ) * inc1
           ic = ib
           ja = 1
@@ -1506,6 +1534,6 @@
                 DO  ijk = 1, lot
                    c(ja+j) = a(ia+i) + a(ib+i)
-                   c(jb+j) = (a(ia+i)-0.5_wp*a(ib+i)) - (sin60*(b(ib+i)))
-                   c(jc+j) = (a(ia+i)-0.5_wp*a(ib+i)) + (sin60*(b(ib+i)))
+                   c(jb+j) = ( a(ia+i) - 0.5_wp * a(ib+i) ) - ( sin60 * ( b(ib+i) ) )
+                   c(jc+j) = ( a(ia+i) - 0.5_wp * a(ib+i) ) + ( sin60 * ( b(ib+i) ) )
                    i = i + inc3
                    j = j + inc4
@@ -1538,14 +1566,18 @@
                       !OCL NOVREC
                       DO  ijk = 1, lot
-                         c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
-                         d(ja+j) = b(ia+i) + (b(ib+i)-b(ic+i))
-                         c(jb+j) = c1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ b(ic+i)))) &
-                                   - s1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- a(ic+i))))
-                         d(jb+j) = s1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ b(ic+i)))) &
-                                   + c1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- a(ic+i))))
-                         c(jc+j) = c2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ b(ic+i)))) &
-                                   - s2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- a(ic+i))))
-                         d(jc+j) = s2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ b(ic+i)))) &
-                                   + c2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- a(ic+i))))
+                         c(ja+j) = a(ia+i) + ( a(ib+i) + a(ic+i) )
+                         d(ja+j) = b(ia+i) + ( b(ib+i) - b(ic+i) )
+                         c(jb+j) = c1 * ( ( a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) ) ) - ( sin60   &
+                                   * ( b(ib+i) + b(ic+i) ) ) ) - s1 * ( ( b(ia+i) - 0.5_wp         &
+                                   * ( b(ib+i) - b(ic+i) ) ) + ( sin60 * ( a(ib+i) - a(ic+i) ) ) )
+                         d(jb+j) = s1 * ( ( a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) ) ) - ( sin60   &
+                                   * ( b(ib+i) + b(ic+i) ) ) ) + c1 * ( ( b(ia+i) - 0.5_wp         &
+                                   * ( b(ib+i) - b(ic+i) ) ) + ( sin60 * ( a(ib+i) - a(ic+i) ) ) )
+                         c(jc+j) = c2 * ( ( a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) ) ) + ( sin60   &
+                                   * ( b(ib+i) + b(ic+i) ) ) ) - s2 * ( ( b(ia+i) - 0.5_wp         &
+                                   * ( b(ib+i) - b(ic+i) ) ) - ( sin60 * ( a(ib+i) - a(ic+i) ) ) )
+                         d(jc+j) = s2 * ( ( a(ia+i) - 0.5_wp * ( a(ib+i) + a(ic+i) ) ) + ( sin60   &
+                                   * ( b(ib+i) + b(ic+i) ) ) ) + c2 * ( ( b(ia+i) - 0.5_wp         &
+                                   * ( b(ib+i) - b(ic+i) ) ) - ( sin60 * ( a(ib+i) - a(ic+i) ) ) )
                          i = i + inc3
                          j = j + inc4
@@ -1573,6 +1605,6 @@
                 DO  ijk = 1, lot
                    c(ja+j) = a(ia+i) + a(ib+i)
-                   c(jb+j) = (0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
-                   c(jc+j) = -(0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
+                   c(jb+j) = ( 0.5_wp * a(ia+i) - a(ib+i) ) - ( sin60 * b(ia+i) )
+                   c(jc+j) = - ( 0.5_wp * a(ia+i) - a(ib+i) ) - ( sin60 * b(ia+i) )
                    i = i + inc3
                    j = j + inc4
@@ -1592,7 +1624,7 @@
                 !OCL NOVREC
                 DO  ijk = 1, lot
-                   c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i))
-                   c(jb+j) = (2.0_wp*a(ia+i)-a(ib+i)) - (ssin60*b(ib+i))
-                   c(jc+j) = (2.0_wp*a(ia+i)-a(ib+i)) + (ssin60*b(ib+i))
+                   c(ja+j) = 2.0_wp * ( a(ia+i) + a(ib+i) )
+                   c(jb+j) = ( 2.0_wp * a(ia+i) - a(ib+i) ) - ( ssin60 * b(ib+i) )
+                   c(jc+j) = ( 2.0_wp * a(ia+i) - a(ib+i) ) + ( ssin60 * b(ib+i) )
                    i = i + inc3
                    j = j + inc4
@@ -1609,6 +1641,6 @@
 
           ia = 1
-          ib = ia + (2*m-la) * inc1
-          ic = ib + 2*m*inc1
+          ib = ia + ( 2 * m - la ) * inc1
+          ic = ib + 2 * m * inc1
           id = ib
           ja = 1
@@ -1626,8 +1658,8 @@
                 !OCL NOVREC
                 DO  ijk = 1, lot
-                   c(ja+j) = (a(ia+i)+a(ic+i)) + a(ib+i)
-                   c(jb+j) = (a(ia+i)-a(ic+i)) - b(ib+i)
-                   c(jc+j) = (a(ia+i)+a(ic+i)) - a(ib+i)
-                   c(jd+j) = (a(ia+i)-a(ic+i)) + b(ib+i)
+                   c(ja+j) = ( a(ia+i) + a(ic+i) ) + a(ib+i)
+                   c(jb+j) = ( a(ia+i) - a(ic+i) ) - b(ib+i)
+                   c(jc+j) = ( a(ia+i) + a(ic+i) ) - a(ib+i)
+                   c(jd+j) = ( a(ia+i) - a(ic+i) ) + b(ib+i)
                    i = i + inc3
                    j = j + inc4
@@ -1636,11 +1668,11 @@
                 jbase = jbase + inc2
              ENDDO
-             ia = ia + iink
-             iink = 2 * iink
-             ib = ib + iink
-             ic = ic - iink
-             id = id - iink
+             ia    = ia + iink
+             iink  = 2 * iink
+             ib    = ib + iink
+             ic    = ic - iink
+             id    = id - iink
              jbase = jbase + jump
-             jump = 2 * jump + jink
+             jump  = 2 * jump + jink
 
              IF ( ib /= ic )  THEN
@@ -1664,18 +1696,18 @@
                       !OCL NOVREC
                       DO  ijk = 1, lot
-                         c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))
-                         d(ja+j) = (b(ia+i)-b(ic+i)) + (b(ib+i)-b(id+i))
-                         c(jc+j) = c2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) - s2*((b(ia+i)-b(ic+i))&
-                                   -(b(ib+i)-b(id+i)))
-                         d(jc+j) = s2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) + c2*((b(ia+i)-b(ic+i))&
-                                   -(b(ib+i)-b(id+i)))
-                         c(jb+j) = c1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) - s1*((b(ia+i)+b(ic+i))&
-                                   +(a(ib+i)-a(id+i)))
-                         d(jb+j) = s1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) + c1*((b(ia+i)+b(ic+i))&
-                                   +(a(ib+i)-a(id+i)))
-                         c(jd+j) = c3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) - s3*((b(ia+i)+b(ic+i))&
-                                   -(a(ib+i)-a(id+i)))
-                         d(jd+j) = s3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) + c3*((b(ia+i)+b(ic+i))&
-                                   -(a(ib+i)-a(id+i)))
+                         c(ja+j) = ( a(ia+i) + a(ic+i) ) + ( a(ib+i) + a(id+i) )
+                         d(ja+j) = ( b(ia+i) - b(ic+i) ) + ( b(ib+i) - b(id+i) )
+                         c(jc+j) = c2 * ( ( a(ia+i) + a(ic+i) ) - ( a(ib+i) + a(id+i) ) ) - s2     &
+                                   * ( ( b(ia+i) - b(ic+i) ) - ( b(ib+i) - b(id+i) ) )
+                         d(jc+j) = s2 * ( ( a(ia+i) + a(ic+i) ) - ( a(ib+i) + a(id+i) ) ) + c2     &
+                                   * ( ( b(ia+i) - b(ic+i) ) - ( b(ib+i) - b(id+i) ) )
+                         c(jb+j) = c1 * ( ( a(ia+i) - a(ic+i) ) - ( b(ib+i) + b(id+i) ) ) - s1     &
+                                   * ( ( b(ia+i) + b(ic+i) ) + ( a(ib+i) - a(id+i) ) )
+                         d(jb+j) = s1 * ( ( a(ia+i) - a(ic+i) ) - ( b (ib+i) + b(id+i) ) ) + c1    &
+                                   * ( ( b(ia+i) + b(ic+i) ) + ( a(ib+i) - a(id+i) ) )
+                         c(jd+j) = c3 * ( ( a(ia+i) - a(ic+i) ) + ( b(ib+i) + b(id+i) ) ) - s3     &
+                                   * ( ( b(ia+i) + b(ic+i) ) - ( a(ib+i) - a(id+i) ) )
+                         d(jd+j) = s3 * ( ( a(ia+i) - a(ic+i) ) + ( b(ib+i) + b(id+i) ) ) + c3     &
+                                   * ( ( b(ia+i) + b(ic+i) ) - ( a(ib+i) - a(id+i) ) )
                          i = i + inc3
                          j = j + inc4
@@ -1705,7 +1737,7 @@
                 DO  ijk = 1, lot
                    c(ja+j) = a(ia+i) + a(ib+i)
-                   c(jb+j) = sin45*((a(ia+i)-a(ib+i))-(b(ia+i)+b(ib+i)))
+                   c(jb+j) = sin45 * ( ( a(ia+i) - a(ib+i) ) - ( b(ia+i) + b(ib+i) ) )
                    c(jc+j) = b(ib+i) - b(ia+i)
-                   c(jd+j) = -sin45*((a(ia+i)-a(ib+i))+(b(ia+i)+b(ib+i)))
+                   c(jd+j) = - sin45 * ( ( a(ia+i) - a(ib+i) ) + ( b(ia+i) + b(ib+i) ) )
                    i = i + inc3
                    j = j + inc4
@@ -1724,8 +1756,8 @@
                 !OCL NOVREC
                 DO  ijk = 1, lot
-                   c(ja+j) = 2.0_wp*((a(ia+i)+a(ic+i))+a(ib+i))
-                   c(jb+j) = 2.0_wp*((a(ia+i)-a(ic+i))-b(ib+i))
-                   c(jc+j) = 2.0_wp*((a(ia+i)+a(ic+i))-a(ib+i))
-                   c(jd+j) = 2.0_wp*((a(ia+i)-a(ic+i))+b(ib+i))
+                   c(ja+j) = 2.0_wp * ( ( a(ia+i) + a(ic+i) ) + a(ib+i) )
+                   c(jb+j) = 2.0_wp * ( ( a(ia+i) - a(ic+i) ) - b(ib+i) )
+                   c(jc+j) = 2.0_wp * ( ( a(ia+i) + a(ic+i) ) - a(ib+i) )
+                   c(jd+j) = 2.0_wp * ( ( a(ia+i) - a(ic+i) ) + b(ib+i) )
                    i = i + inc3
                    j = j + inc4
@@ -1742,6 +1774,6 @@
 
           ia = 1
-          ib = ia + (2*m-la) * inc1
-          ic = ib + 2*m*inc1
+          ib = ia + ( 2 * m - la ) * inc1
+          ic = ib + 2 * m * inc1
           id = ic
           ie = ib
@@ -1761,13 +1793,13 @@
                 !OCL NOVREC
                 DO  ijk = 1, lot
-                   c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
-                   c(jb+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) -        &
-                             (sin72*b(ib+i)+sin36*b(ic+i))
-                   c(jc+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) -        &
-                             (sin36*b(ib+i)-sin72*b(ic+i))
-                   c(jd+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) +        &
-                             (sin36*b(ib+i)-sin72*b(ic+i))
-                   c(je+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) +        &
-                             (sin72*b(ib+i)+sin36*b(ic+i))
+                   c(ja+j) = a(ia+i) + ( a(ib+i) + a(ic+i) )
+                   c(jb+j) = ( ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) + qrt5 *              &
+                             ( a(ib+i) - a(ic+i) ) ) - ( sin72 * b(ib+i) + sin36 * b(ic+i) )
+                   c(jc+j) = ( ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) - qrt5 *              &
+                             ( a(ib+i) - a(ic+i) ) ) - ( sin36 * b(ib+i) - sin72 * b(ic+i) )
+                   c(jd+j) = ( ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) - qrt5 *              &
+                             ( a(ib+i) - a(ic+i) ) ) + ( sin36 * b(ib+i) - sin72 * b(ic+i) )
+                   c(je+j) = ( ( a(ia+i) - 0.25_wp * ( a(ib+i) + a(ic+i) ) ) + qrt5 *              &
+                             ( a(ib+i) - a(ic+i) ) ) + ( sin72 * b(ib+i) + sin36 * b(ic+i) )
                    i = i + inc3
                    j = j + inc4
@@ -1776,12 +1808,12 @@
                 jbase = jbase + inc2
              ENDDO
-             ia = ia + iink
-             iink = 2 * iink
-             ib = ib + iink
-             ic = ic + iink
-             id = id - iink
-             ie = ie - iink
+             ia    = ia + iink
+             iink  = 2 * iink
+             ib    = ib + iink
+             ic    = ic + iink
+             id    = id - iink
+             ie    = ie - iink
              jbase = jbase + jump
-             jump = 2 * jump + jink
+             jump  = 2 * jump + jink
 
              IF ( ib /= id )  THEN
@@ -1808,27 +1840,30 @@
                       !OCL NOVREC
                       DO  ijk = 1, lot
-                         a10(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) +      &
-                                    qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
-                         a20(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) -      &
-                                    qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
-                         b10(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) +      &
-                                    qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
-                         b20(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) -      &
-                                    qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
-                         a11(ijk) = sin72*(b(ib+i)+b(ie+i)) + sin36*(b(ic+i)+b(id+i))
-                         a21(ijk) = sin36*(b(ib+i)+b(ie+i)) - sin72*(b(ic+i)+b(id+i))
-                         b11(ijk) = sin72*(a(ib+i)-a(ie+i)) + sin36*(a(ic+i)-a(id+i))
-                         b21(ijk) = sin36*(a(ib+i)-a(ie+i)) - sin72*(a(ic+i)-a(id+i))
-                         c(ja+j) = a(ia+i) + ((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))
-                         d(ja+j) = b(ia+i) + ((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))
-                         c(jb+j) = c1*(a10(ijk)-a11(ijk)) - s1*(b10(ijk)+b11(ijk))
-                         d(jb+j) = s1*(a10(ijk)-a11(ijk)) + c1*(b10(ijk)+b11(ijk))
-                         c(je+j) = c4*(a10(ijk)+a11(ijk)) - s4*(b10(ijk)-b11(ijk))
-                         d(je+j) = s4*(a10(ijk)+a11(ijk)) + c4*(b10(ijk)-b11(ijk))
-                         c(jc+j) = c2*(a20(ijk)-a21(ijk)) - s2*(b20(ijk)+b21(ijk))
-                         d(jc+j) = s2*(a20(ijk)-a21(ijk)) + c2*(b20(ijk)+b21(ijk))
-                         c(jd+j) = c3*(a20(ijk)+a21(ijk)) - s3*(b20(ijk)-b21(ijk))
-                         d(jd+j) = s3*(a20(ijk)+a21(ijk)) + c3*(b20(ijk)-b21(ijk))
-
+                         a10(ijk) = ( a(ia+i) - 0.25_wp * ( ( a(ib+i) + a(ie+i) ) + ( a(ic+i)      &
+                                    + a(id+i) ) ) ) + qrt5 * ( ( a(ib+i) + a(ie+i) ) - ( a(ic+i)   &
+                                    + a(id+i) ) )
+                         a20(ijk) = ( a(ia+i) - 0.25_wp * ( ( a(ib+i) + a(ie+i) ) + ( a(ic+i)      &
+                                    + a(id+i) ) ) ) - qrt5 * ( ( a(ib+i) + a(ie+i) ) - ( a(ic+i)   &
+                                    + a(id+i) ) )
+                         b10(ijk) = ( b(ia+i) - 0.25_wp * ( ( b(ib+i) - b(ie+i) ) + ( b(ic+i)      &
+                                    - b(id+i) ) ) ) + qrt5 * ( ( b(ib+i) - b(ie+i) ) - ( b(ic+i)   &
+                                    - b(id+i) ) )
+                         b20(ijk) = ( b(ia+i) - 0.25_wp * ( ( b(ib+i) - b(ie+i) ) + ( b(ic+i)      &
+                                    - b(id+i) ) ) ) - qrt5 * ( ( b(ib+i) - b(ie+i) ) - ( b(ic+i)   &
+                                    - b(id+i) ) )
+                         a11(ijk) = sin72 * ( b(ib+i) + b(ie+i) ) + sin36 * ( b(ic+i) + b(id+i) )
+                         a21(ijk) = sin36 * ( b(ib+i) + b(ie+i) ) - sin72 * ( b(ic+i) + b(id+i) )
+                         b11(ijk) = sin72 * ( a(ib+i) - a(ie+i) ) + sin36 * ( a(ic+i) - a(id+i) )
+                         b21(ijk) = sin36 * ( a(ib+i) - a(ie+i) ) - sin72 * ( a(ic+i) - a(id+i) )
+                         c(ja+j)  = a(ia+i) + ( ( a(ib+i) + a(ie+i) ) + ( a(ic+i) + a(id+i) ) )
+                         d(ja+j)  = b(ia+i) + ( ( b(ib+i) - b(ie+i) ) + ( b(ic+i) - b(id+i) ) )
+                         c(jb+j)  = c1 * ( a10(ijk) - a11(ijk) ) - s1 * ( b10(ijk) + b11(ijk) )
+                         d(jb+j)  = s1 * ( a10(ijk) - a11(ijk) ) + c1 * ( b10(ijk) + b11(ijk) )
+                         c(je+j)  = c4 * ( a10(ijk) + a11(ijk) ) - s4 * ( b10(ijk) - b11(ijk) )
+                         d(je+j)  = s4 * ( a10(ijk) + a11(ijk) ) + c4 * ( b10(ijk) - b11(ijk) )
+                         c(jc+j)  = c2 * ( a20(ijk) - a21(ijk) ) - s2 * ( b20(ijk) + b21(ijk) )
+                         d(jc+j)  = s2 * ( a20(ijk) - a21(ijk) ) + c2 * ( b20(ijk) + b21(ijk) )
+                         c(jd+j)  = c3 * ( a20(ijk) + a21(ijk) ) - s3 * ( b20(ijk) - b21(ijk) )
+                         d(jd+j)  = s3 * ( a20(ijk) + a21(ijk) ) + c3 * ( b20(ijk) - b21(ijk) )
                          i = i + inc3
                          j = j + inc4
@@ -1857,13 +1892,13 @@
                 !OCL NOVREC
                 DO  ijk = 1, lot
-                   c(ja+j) = (a(ia+i)+a(ib+i)) + a(ic+i)
-                   c(jb+j) = (qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) -        &
-                             (sin36*b(ia+i)+sin72*b(ib+i))
-                   c(je+j) = -(qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) -       &
-                             (sin36*b(ia+i)+sin72*b(ib+i))
-                   c(jc+j) = (qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) -        &
-                             (sin72*b(ia+i)-sin36*b(ib+i))
-                   c(jd+j) = -(qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) -       &
-                             (sin72*b(ia+i)-sin36*b(ib+i))
+                   c(ja+j) = ( a(ia+i) + a(ib+i) ) + a(ic+i)
+                   c(jb+j) = ( qrt5 * ( a(ia+i) - a(ib+i) ) + ( 0.25_wp * ( a(ia+i) + a(ib+i) )    &
+                             - a(ic+i) ) ) - ( sin36 * b(ia+i) + sin72 * b(ib+i) )
+                   c(je+j) = - ( qrt5 * ( a(ia+i) - a(ib+i) ) + ( 0.25_wp * ( a(ia+i) + a(ib+i) )  &
+                             - a(ic+i) ) ) - ( sin36 * b(ia+i) + sin72 * b(ib+i) )
+                   c(jc+j) = ( qrt5 * ( a(ia+i) - a(ib+i) ) - ( 0.25_wp * ( a(ia+i) + a(ib+i) )    &
+                             - a(ic+i) ) ) - ( sin72 * b(ia+i) - sin36 * b(ib+i) )
+                   c(jd+j) = - ( qrt5 * ( a(ia+i) - a(ib+i) ) - ( 0.25_wp * ( a(ia+i) + a(ib+i) )  &
+                             - a(ic+i) ) ) - ( sin72 * b(ia+i) - sin36 * b(ib+i) )
                    i = i + inc3
                    j = j + inc4
@@ -1885,13 +1920,13 @@
                 !OCL NOVREC
                 DO  ijk = 1, lot
-                   c(ja+j) = 2.0_wp*(a(ia+i)+(a(ib+i)+a(ic+i)))
-                   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)))  &
-                             - (ssin72*b(ib+i)+ssin36*b(ic+i))
-                   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)))  &
-                             - (ssin36*b(ib+i)-ssin72*b(ic+i))
-                   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)))  &
-                             + (ssin36*b(ib+i)-ssin72*b(ic+i))
-                   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)))  &
-                             + (ssin72*b(ib+i)+ssin36*b(ic+i))
+                   c(ja+j) = 2.0_wp * ( a(ia+i) + ( a(ib+i) + a(ic+i) ) )
+                   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) ) ) - ( ssin72 * b(ib+i) + ssin36 * b(ic+i) )
+                   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) ) ) - ( ssin36 * b(ib+i) - ssin72 * b(ic+i) )
+                   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) ) ) + ( ssin36 * b(ib+i) - ssin72 * b(ic+i) )
+                   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) ) ) + (ssin72 * b(ib+i) + ssin36 * b(ic+i) )
                    i = i + inc3
                    j = j + inc4
@@ -1908,7 +1943,7 @@
 
           ia = 1
-          ib = ia + (2*m-la) * inc1
-          ic = ib + 2*m*inc1
-          id = ic + 2*m*inc1
+          ib = ia + ( 2 * m - la ) * inc1
+          ic = ib + 2 * m * inc1
+          id = ic + 2 * m * inc1
           ie = ic
           if = ib
@@ -1929,10 +1964,14 @@
                 !OCL NOVREC
                 DO  ijk = 1, lot
-                   c(ja+j) = (a(ia+i)+a(id+i)) + (a(ib+i)+a(ic+i))
-                   c(jd+j) = (a(ia+i)-a(id+i)) - (a(ib+i)-a(ic+i))
-                   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)))
-                   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)))
-                   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)))
-                   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)))
+                   c(ja+j) = ( a (ia+i) + a(id+i) ) +  ( a(ib+i) + a(ic+i) )
+                   c(jd+j) = ( a (ia+i) - a(id+i) ) -  ( a(ib+i) - a(ic+i) )
+                   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) ) )
+                   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) )
+                   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) )
+                   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) )
                    i = i + inc3
                    j = j + inc4
@@ -1941,13 +1980,13 @@
                 jbase = jbase + inc2
              ENDDO
-             ia = ia + iink
-             iink = 2 * iink
-             ib = ib + iink
-             ic = ic + iink
-             id = id - iink
-             ie = ie - iink
-             if = if - iink
+             ia    = ia + iink
+             iink  = 2 * iink
+             ib    = ib + iink
+             ic    = ic + iink
+             id    = id - iink
+             ie    = ie - iink
+             if    = if - iink
              jbase = jbase + jump
-             jump = 2 * jump + jink
+             jump  = 2 * jump + jink
 
              IF ( ic /= id )  THEN
@@ -1977,31 +2016,33 @@
                       !OCL NOVREC
                       DO  ijk = 1, lot
-                         a11(ijk) = (a(ie+i)+a(ib+i)) + (a(ic+i)+a(if+i))
-                         a20(ijk) = (a(ia+i)+a(id+i)) - 0.5_wp*a11(ijk)
-                         a21(ijk) = sin60*((a(ie+i)+a(ib+i))-(a(ic+i)+a(if+i)))
-                         b11(ijk) = (b(ib+i)-b(ie+i)) + (b(ic+i)-b(if+i))
-                         b20(ijk) = (b(ia+i)-b(id+i)) - 0.5_wp*b11(ijk)
-                         b21(ijk) = sin60*((b(ib+i)-b(ie+i))-(b(ic+i)-b(if+i)))
-
-                         c(ja+j) = (a(ia+i)+a(id+i)) + a11(ijk)
-                         d(ja+j) = (b(ia+i)-b(id+i)) + b11(ijk)
-                         c(jc+j) = c2*(a20(ijk)-b21(ijk)) - s2*(b20(ijk)+a21(ijk))
-                         d(jc+j) = s2*(a20(ijk)-b21(ijk)) + c2*(b20(ijk)+a21(ijk))
-                         c(je+j) = c4*(a20(ijk)+b21(ijk)) - s4*(b20(ijk)-a21(ijk))
-                         d(je+j) = s4*(a20(ijk)+b21(ijk)) + c4*(b20(ijk)-a21(ijk))
-
-                         a11(ijk) = (a(ie+i)-a(ib+i)) + (a(ic+i)-a(if+i))
-                         b11(ijk) = (b(ie+i)+b(ib+i)) - (b(ic+i)+b(if+i))
-                         a20(ijk) = (a(ia+i)-a(id+i)) - 0.5_wp*a11(ijk)
-                         a21(ijk) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
-                         b20(ijk) = (b(ia+i)+b(id+i)) + 0.5_wp*b11(ijk)
-                         b21(ijk) = sin60*((b(ie+i)+b(ib+i))+(b(ic+i)+b(if+i)))
-
-                         c(jd+j) = c3*((a(ia+i)-a(id+i))+a11(ijk)) - s3*((b(ia+i)+b(id+i))-b11(ijk))
-                         d(jd+j) = s3*((a(ia+i)-a(id+i))+a11(ijk)) + c3*((b(ia+i)+b(id+i))-b11(ijk))
-                         c(jb+j) = c1*(a20(ijk)-b21(ijk)) - s1*(b20(ijk)-a21(ijk))
-                         d(jb+j) = s1*(a20(ijk)-b21(ijk)) + c1*(b20(ijk)-a21(ijk))
-                         c(jf+j) = c5*(a20(ijk)+b21(ijk)) - s5*(b20(ijk)+a21(ijk))
-                         d(jf+j) = s5*(a20(ijk)+b21(ijk)) + c5*(b20(ijk)+a21(ijk))
+                         a11(ijk) = ( a(ie+i) + a(ib+i) ) + ( a(ic+i) + a(if+i) )
+                         a20(ijk) = ( a(ia+i) + a(id+i) ) - 0.5_wp * a11(ijk)
+                         a21(ijk) = sin60 * ( ( a(ie+i) + a(ib+i) ) - ( a(ic+i) + a(if+i) ) )
+                         b11(ijk) = ( b(ib+i) - b(ie+i) ) + ( b(ic+i) - b(if+i) )
+                         b20(ijk) = ( b(ia+i) - b(id+i) ) - 0.5_wp * b11(ijk)
+                         b21(ijk) = sin60 * ( ( b(ib+i) - b(ie+i) ) - ( b(ic+i) - b(if+i) ) )
+
+                         c(ja+j) = ( a(ia+i) + a(id+i) ) + a11(ijk)
+                         d(ja+j) = ( b(ia+i) - b(id+i) ) + b11(ijk)
+                         c(jc+j) = c2 * ( a20(ijk) - b21(ijk) ) - s2 * ( b20(ijk) + a21(ijk) )
+                         d(jc+j) = s2 * ( a20(ijk) - b21(ijk) ) + c2 * ( b20(ijk) + a21(ijk) )
+                         c(je+j) = c4 * ( a20(ijk) + b21(ijk) ) - s4 * ( b20(ijk) - a21(ijk) )
+                         d(je+j) = s4 * ( a20(ijk) + b21(ijk) ) + c4 * ( b20(ijk) - a21(ijk) )
+
+                         a11(ijk) = ( a(ie+i) - a(ib+i) ) + ( a(ic+i) - a(if+i) )
+                         b11(ijk) = ( b(ie+i) + b(ib+i) ) - ( b(ic+i) + b(if+i) )
+                         a20(ijk) = ( a(ia+i) - a(id+i) ) - 0.5_wp * a11(ijk)
+                         a21(ijk) = sin60 * ( ( a(ie+i) - a(ib+i) ) - ( a(ic+i) - a(if+i) ) )
+                         b20(ijk) = ( b(ia+i) + b(id+i) ) + 0.5_wp * b11(ijk)
+                         b21(ijk) = sin60 * ( ( b(ie+i) + b(ib+i) ) + ( b(ic+i) + b(if+i) ) )
+
+                         c(jd+j) = c3 * ( (a(ia+i) - a(id+i) ) + a11(ijk) ) - s3                   &
+                                      * ( (b(ia+i) + b(id+i) ) - b11(ijk) )
+                         d(jd+j) = s3 * ( (a(ia+i) - a(id+i) ) + a11(ijk) ) + c3                   &
+                                      * ( (b(ia+i) + b(id+i) ) - b11(ijk) )
+                         c(jb+j) = c1 * ( a20(ijk) - b21(ijk) ) - s1 * ( b20(ijk) - a21(ijk) )
+                         d(jb+j) = s1 * ( a20(ijk) - b21(ijk) ) + c1 * ( b20(ijk) - a21(ijk) )
+                         c(jf+j) = c5 * ( a20(ijk) + b21(ijk) ) - s5 * ( b20(ijk) + a21(ijk) )
+                         d(jf+j) = s5 * ( a20(ijk) + b21(ijk) ) + c5 * ( b20(ijk) + a21(ijk) )
 
                          i = i + inc3
@@ -2032,10 +2073,14 @@
                 !OCL NOVREC
                 DO  ijk = 1, lot
-                   c(ja+j) = a(ib+i) + (a(ia+i)+a(ic+i))
-                   c(jd+j) = b(ib+i) - (b(ia+i)+b(ic+i))
-                   c(jb+j) = (sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i))
-                   c(jf+j) = -(sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i))
-                   c(jc+j) = sin60*(b(ic+i)-b(ia+i)) + (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i))
-                   c(je+j) = sin60*(b(ic+i)-b(ia+i)) - (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i))
+                   c(ja+j) = a(ib+i) + ( a(ia+i) + a(ic+i) )
+                   c(jd+j) = b(ib+i) - ( b(ia+i) + b(ic+i) )
+                   c(jb+j) = ( sin60 * ( a(ia+i) - a(ic+i) ) ) - ( 0.5_wp * ( b(ia+i) + b(ic+i) )  &
+                             + b(ib+i) )
+                   c(jf+j) = - ( sin60 * ( a(ia+i) - a(ic+i) ) ) - ( 0.5_wp * ( b(ia+i)            &
+                             + b(ic+i) ) + b(ib+i) )
+                   c(jc+j) = sin60 * ( b(ic+i) - b(ia+i) ) + ( 0.5_wp * ( a(ia+i) + a(ic+i) )      &
+                             - a(ib+i) )
+                   c(je+j) = sin60 * ( b(ic+i) - b(ia+i) ) - ( 0.5_wp * ( a(ia+i) + a(ic+i) )      &
+                             - a(ib+i) )
                    i = i + inc3
                    j = j + inc4
@@ -2055,10 +2100,14 @@
                 !OCL NOVREC
                 DO  ijk = 1, lot
-                   c(ja+j) = (2.0_wp*(a(ia+i)+a(id+i))) + (2.0_wp*(a(ib+i)+a(ic+i)))
-                   c(jd+j) = (2.0_wp*(a(ia+i)-a(id+i))) - (2.0_wp*(a(ib+i)-a(ic+i)))
-                   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)))
-                   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)))
-                   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)))
-                   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)))
+                   c(ja+j) = ( 2.0_wp * ( a(ia+i) + a(id+i) ) ) + ( 2.0_wp * ( a(ib+i) + a(ic+i) ) )
+                   c(jd+j) = ( 2.0_wp * ( a(ia+i) - a(id+i) ) ) - ( 2.0_wp * ( a(ib+i) - a(ic+i) ) )
+                   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) ) )
+                   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) ) )
+                   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) ) )
+                   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) ) )
                    i = i + inc3
                    j = j + inc4
@@ -2080,8 +2129,8 @@
 
           ia = 1
-          ib = ia + la*inc1
-          ic = ib + 2*la*inc1
-          id = ic + 2*la*inc1
-          ie = id + 2*la*inc1
+          ib = ia + la * inc1
+          ic = ib + 2 * la * inc1
+          id = ic + 2 * la * inc1
+          ie = id + 2 * la * inc1
           ja = 1
           jb = ja + jink
@@ -2101,12 +2150,16 @@
              !OCL NOVREC
              DO  ijk = 1, lot
-                c(ja+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))+(a(ib+i)+a(id+i)))
-                c(je+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))-(a(ib+i)+a(id+i)))
-                c(jc+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))-(b(ib+i)-b(id+i)))
-                c(jg+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))+(b(ib+i)-b(id+i)))
-                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)))
-                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)))
-                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)))
-                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)))
+                c(ja+j) = 2.0_wp * ( ( ( a(ia+i) + a(ie+i) ) + a(ic+i) ) + ( a(ib+i) + a(id+i) ) )
+                c(je+j) = 2.0_wp * ( ( ( a(ia+i) + a(ie+i) ) + a(ic+i) ) - ( a(ib+i) + a(id+i) ) )
+                c(jc+j) = 2.0_wp * ( ( ( a(ia+i) + a(ie+i) ) - a(ic+i) ) - ( b(ib+i) - b(id+i) ) )
+                c(jg+j) = 2.0_wp * ( ( ( a(ia+i) + a(ie+i) ) - a(ic+i) ) + ( b(ib+i) - b(id+i) ) )
+                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) ) )
+                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) ) )
+                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) ) )
+                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) ) )
                 i = i + inc3
                 j = j + inc4
@@ -2120,12 +2173,11 @@
  END SUBROUTINE rpassm
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Computes factors of n & trigonometric functins required by fft99 & fft991cy
 !> Method: Dimension trigs(n),ifax(1),jfax(10),lfax(7)
-!> subroutine 'set99' - computes factors of n & trigonometric
-!> functins required by fft99 & fft991cy
-!------------------------------------------------------------------------------!
+!> Subroutine 'set99' - computes factors of n & trigonometric functins required by fft99 & fft991cy
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE set99( trigs, ifax, n )
 
@@ -2138,24 +2190,23 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  n !< 
-
-    INTEGER(iwp) ::  ifax(*)  !< 
-    REAL(wp)     ::  trigs(*) !< 
-
-
-    REAL(wp) ::  angle    !< 
-    REAL(wp) ::  del      !< 
-    INTEGER(iwp) ::  i    !< 
-    INTEGER(iwp) ::  ifac !< 
-    INTEGER(iwp) ::  ixxx !< 
-    INTEGER(iwp) ::  k    !< 
-    INTEGER(iwp) ::  l    !< 
-    INTEGER(iwp) ::  nfax !< 
-    INTEGER(iwp) ::  nhl  !< 
-    INTEGER(iwp) ::  nil  !< 
-    INTEGER(iwp) ::  nu   !< 
-
-    INTEGER(iwp) ::  jfax(10) !< 
-    INTEGER(iwp) ::  lfax(7)  !< 
+    INTEGER(iwp) ::  ifax(*)   !<
+
+    INTEGER(iwp) ::  i         !<
+    INTEGER(iwp) ::  ifac      !<
+    INTEGER(iwp) ::  ixxx      !<
+    INTEGER(iwp) ::  jfax(10)  !<
+    INTEGER(iwp) ::  k         !<
+    INTEGER(iwp) ::  l         !<
+    INTEGER(iwp) ::  lfax(7)   !<
+    INTEGER(iwp) ::  n         !<
+    INTEGER(iwp) ::  nfax      !<
+    INTEGER(iwp) ::  nhl       !<
+    INTEGER(iwp) ::  nil       !<
+    INTEGER(iwp) ::  nu        !<
+
+    REAL(wp) ::  angle         !<
+    REAL(wp) ::  del           !<
+
+    REAL(wp) ::  trigs(*)      !<
 
     DATA  lfax / 6, 8, 5, 4, 3, 2, 1 /
@@ -2164,9 +2215,9 @@
     ixxx = 1
 
-    del = 4.0_wp * ASIN( 1.0_wp ) / REAL( n, KIND=wp )
+    del = 4.0_wp * ASIN( 1.0_wp ) / REAL( n, KIND = wp )
     nil = 0
-    nhl = (n/2) - 1
+    nhl = ( n / 2 ) - 1
     DO  k = nil, nhl
-       angle = REAL( k, KIND=wp ) * del
+       angle = REAL( k, KIND = wp ) * del
        trigs(2*k+1) = COS( angle )
        trigs(2*k+2) = SIN( angle )
@@ -2198,6 +2249,5 @@
        ifac = lfax(l)
        IF (ifac < 2 )  THEN
-          message_string = 'number of gridpoints along x or/and y ' //                             &
-                           'contain illegal  factors' //                                           &
+          message_string = 'number of gridpoints along x or/and y contain illegal  factors' //     &
                            '&only factors 2, 3, 5 are allowed'
           CALL message( 'temperton_fft', 'PA0311', 1, 2, 0, 6, 0 )
@@ -2224,33 +2274,33 @@
     IMPLICIT NONE
 
-    REAL(wp),DIMENSION(:,:)     ::  a    !<
-    REAL(wp),DIMENSION(:)       ::  trigs !<
-    REAL(wp),DIMENSION(:,:)     ::  work  !<
+
+    INTEGER(iwp) ::  i       !<
+    INTEGER(iwp) ::  ia      !<
+    INTEGER(iwp) ::  ibase   !<
+    INTEGER(iwp) ::  ierr    !<
+    INTEGER(iwp) ::  ifac    !<
+    INTEGER(iwp) ::  igo     !<
+    INTEGER(iwp) ::  ii      !<
+    INTEGER(iwp) ::  inc     !<
+    INTEGER(iwp) ::  isign   !<
+    INTEGER(iwp) ::  istart  !<
+    INTEGER(iwp) ::  ix      !<
+    INTEGER(iwp) ::  iz      !<
+    INTEGER(iwp) ::  j       !<
+    INTEGER(iwp) ::  jbase   !<
+    INTEGER(iwp) ::  jump    !<
+    INTEGER(iwp) ::  k       !<
+    INTEGER(iwp) ::  la      !<
+    INTEGER(iwp) ::  lot     !<
+    INTEGER(iwp) ::  n       !<
+    INTEGER(iwp) ::  nfax    !<
+    INTEGER(iwp) ::  nvex    !<
+    INTEGER(iwp) ::  nx      !<
+
     INTEGER(iwp),DIMENSION(:),INTENT(IN) ::  ifax  !<
 
-    INTEGER(iwp) ::  inc   !<
-    INTEGER(iwp) ::  isign !<
-    INTEGER(iwp) ::  jump  !<
-    INTEGER(iwp) ::  lot   !<
-    INTEGER(iwp) ::  n     !<
-
-
-    INTEGER(iwp) ::  i      !<
-    INTEGER(iwp) ::  ia     !<
-    INTEGER(iwp) ::  ibase  !<
-    INTEGER(iwp) ::  ierr   !<
-    INTEGER(iwp) ::  ifac   !<
-    INTEGER(iwp) ::  igo    !<
-    INTEGER(iwp) ::  ii     !<
-    INTEGER(iwp) ::  istart !<
-    INTEGER(iwp) ::  ix     !<
-    INTEGER(iwp) ::  iz     !<
-    INTEGER(iwp) ::  j      !<
-    INTEGER(iwp) ::  jbase  !<
-    INTEGER(iwp) ::  k      !<
-    INTEGER(iwp) ::  la     !<
-    INTEGER(iwp) ::  nfax   !<
-    INTEGER(iwp) ::  nvex   !<
-    INTEGER(iwp) ::  nx     !<
+    REAL(wp),DIMENSION(:,:) ::  a      !<
+    REAL(wp),DIMENSION(:)   ::  trigs  !<
+    REAL(wp),DIMENSION(:,:) ::  work   !<
 
 
@@ -2262,5 +2312,5 @@
     nfax = ifax(1)
     nx   = n + 1
-    IF ( MOD(n,2) == 1 )  nx = n
+    IF ( MOD( n, 2 ) == 1 )  nx = n
     nvex = 1
 
@@ -2273,5 +2323,5 @@
        i = istart
        a(:,i+inc) = 0.5_wp * a(:,i)
-       IF ( MOD(n,2) /= 1 )  THEN
+       IF ( MOD( n, 2 ) /= 1 )  THEN
           i = istart + n * inc
           a(:,i) = 0.5_wp * a(:,i)
@@ -2313,5 +2363,5 @@
 !
 !--    If necessary, copy results back to a
-       IF ( MOD(nfax,2) /= 0 )  THEN
+       IF ( MOD( nfax, 2 ) /= 0 )  THEN
           ibase = 1
           jbase = ia
@@ -2326,5 +2376,5 @@
 !
 !--    Fill in zeros at end
-       ix = istart + n*inc
+       ix = istart + n * inc
        a(:,ix) = 0.0_wp
        a(:,ix+inc) = 0.0_wp
@@ -2370,5 +2420,5 @@
 !
 !--    If necessary, copy results back to a
-       IF ( MOD(nfax,2) /= 0 )  THEN
+       IF ( MOD( nfax, 2 ) /= 0 )  THEN
           ibase = 1
           jbase = ia
@@ -2387,6 +2437,6 @@
        a(:,ix+inc) = 0.0_wp
 
-       IF ( MOD(n,2) /= 1 )  THEN
-          iz = istart + (n+1) * inc
+       IF ( MOD( n, 2 ) /= 1 )  THEN
+          iz = istart + ( n + 1 ) * inc
           a(:,iz) = 0.0_wp
        ENDIF
@@ -2396,9 +2446,8 @@
  END SUBROUTINE fft991cy_vec
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> Performs one pass through data as part of
-!> multiple real fft (fourier analysis) routine.
+!> Performs one pass through data as part of multiple real fft (fourier analysis) routine.
 !>
 !> Method:
@@ -2422,5 +2471,5 @@
 !>         2 - ifac not catered for
 !>         3 - ifac only catered for if la=n/ifac
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE qpassm_vec( a, b, c, d, trigs, inc1, inc2, n, ifac, la, ierr )
 
@@ -2429,95 +2478,98 @@
     IMPLICIT NONE
 
-    INTEGER(iwp),INTENT(IN)  ::  ifac !<
-    INTEGER(iwp),INTENT(IN)  ::  inc1 !<
-    INTEGER(iwp),INTENT(IN)  ::  inc2 !<
-    INTEGER(iwp),INTENT(IN)  ::  la   !<
-    INTEGER(iwp),INTENT(IN)  ::  n    !<
-    INTEGER(iwp),INTENT(OUT) ::  ierr !<
+    INTEGER(iwp),INTENT(IN)  ::  ifac  !<
+    INTEGER(iwp),INTENT(IN)  ::  inc1  !<
+    INTEGER(iwp),INTENT(IN)  ::  inc2  !<
+    INTEGER(iwp),INTENT(IN)  ::  la    !<
+    INTEGER(iwp),INTENT(IN)  ::  n     !<
+
+    INTEGER(iwp),INTENT(OUT) ::  ierr  !<
+
+    INTEGER(iwp) ::  i      !<
+    INTEGER(iwp) ::  ia     !<
+    INTEGER(iwp) ::  ib     !<
+    INTEGER(iwp) ::  ibase  !<
+    INTEGER(iwp) ::  ic     !<
+    INTEGER(iwp) ::  id     !<
+    INTEGER(iwp) ::  ie     !<
+    INTEGER(iwp) ::  if     !<
+    INTEGER(iwp) ::  ig     !<
+    INTEGER(iwp) ::  igo    !<
+    INTEGER(iwp) ::  ih     !<
+    INTEGER(iwp) ::  iink   !<
+    INTEGER(iwp) ::  ijump  !<
+    INTEGER(iwp) ::  ipl    !<  loop index parallel loop
+    INTEGER(iwp) ::  j      !<
+    INTEGER(iwp) ::  ja     !<
+    INTEGER(iwp) ::  jb     !<
+    INTEGER(iwp) ::  jbase  !<
+    INTEGER(iwp) ::  jc     !<
+    INTEGER(iwp) ::  jd     !<
+    INTEGER(iwp) ::  je     !<
+    INTEGER(iwp) ::  jf     !<
+    INTEGER(iwp) ::  jink   !<
+    INTEGER(iwp) ::  k      !<
+    INTEGER(iwp) ::  kb     !<
+    INTEGER(iwp) ::  kc     !<
+    INTEGER(iwp) ::  kd     !<
+    INTEGER(iwp) ::  ke     !<
+    INTEGER(iwp) ::  kf     !<
+    INTEGER(iwp) ::  kstop  !<
+    INTEGER(iwp) ::  l      !<
+    INTEGER(iwp) ::  m      !<
 
 !
 !-- Arrays are dimensioned with n
-    REAL(wp),DIMENSION(:,:) ::  a     !<
-    REAL(wp),DIMENSION(:,:) ::  b     !<
-    REAL(wp),DIMENSION(:,:) ::  c     !<
-    REAL(wp),DIMENSION(:,:) ::  d     !<
-    REAL(wp),DIMENSION(:),INTENT(IN) ::  trigs !<
-
-    REAL(wp) ::  a0     !<
-    REAL(wp) ::  a1     !<
-    REAL(wp) ::  a10    !<
-    REAL(wp) ::  a11    !<
-    REAL(wp) ::  a2     !<
-    REAL(wp) ::  a20    !<
-    REAL(wp) ::  a21    !<
-    REAL(wp) ::  a3     !<
-    REAL(wp) ::  a4     !<
-    REAL(wp) ::  a5     !<
-    REAL(wp) ::  a6     !<
-    REAL(wp) ::  b0     !<
-    REAL(wp) ::  b1     !<
-    REAL(wp) ::  b10    !<
-    REAL(wp) ::  b11    !<
-    REAL(wp) ::  b2     !<
-    REAL(wp) ::  b20    !<
-    REAL(wp) ::  b21    !<
-    REAL(wp) ::  b3     !<
-    REAL(wp) ::  b4     !<
-    REAL(wp) ::  b5     !<
-    REAL(wp) ::  b6     !<
-    REAL(wp) ::  c1     !<
-    REAL(wp) ::  c2     !<
-    REAL(wp) ::  c3     !<
-    REAL(wp) ::  c4     !<
-    REAL(wp) ::  c5     !<
-    REAL(wp) ::  qrt5   !<
-    REAL(wp) ::  s1     !<
-    REAL(wp) ::  s2     !<
-    REAL(wp) ::  s3     !<
-    REAL(wp) ::  s4     !<
-    REAL(wp) ::  s5     !<
-    REAL(wp) ::  sin36  !<
-    REAL(wp) ::  sin45  !<
-    REAL(wp) ::  sin60  !<
-    REAL(wp) ::  sin72  !<
-    REAL(wp) ::  z      !<
-    REAL(wp) ::  zqrt5  !<
-    REAL(wp) ::  zsin36 !<
-    REAL(wp) ::  zsin45 !<
-    REAL(wp) ::  zsin60 !<
-    REAL(wp) ::  zsin72 !<
-
-    INTEGER(iwp) ::  i     !<
-    INTEGER(iwp) ::  ia    !<
-    INTEGER(iwp) ::  ib    !<
-    INTEGER(iwp) ::  ibase !<
-    INTEGER(iwp) ::  ic    !<
-    INTEGER(iwp) ::  id    !<
-    INTEGER(iwp) ::  ie    !<
-    INTEGER(iwp) ::  if    !<
-    INTEGER(iwp) ::  ig    !<
-    INTEGER(iwp) ::  igo   !<
-    INTEGER(iwp) ::  ih    !<
-    INTEGER(iwp) ::  iink  !<
-    INTEGER(iwp) ::  ijump !<
-    INTEGER(iwp) ::  ipl   !<  loop index parallel loop
-    INTEGER(iwp) ::  j     !<
-    INTEGER(iwp) ::  ja    !<
-    INTEGER(iwp) ::  jb    !<
-    INTEGER(iwp) ::  jbase !<
-    INTEGER(iwp) ::  jc    !<
-    INTEGER(iwp) ::  jd    !<
-    INTEGER(iwp) ::  je    !<
-    INTEGER(iwp) ::  jf    !<
-    INTEGER(iwp) ::  jink  !<
-    INTEGER(iwp) ::  k     !<
-    INTEGER(iwp) ::  kb    !<
-    INTEGER(iwp) ::  kc    !<
-    INTEGER(iwp) ::  kd    !<
-    INTEGER(iwp) ::  ke    !<
-    INTEGER(iwp) ::  kf    !<
-    INTEGER(iwp) ::  kstop !<
-    INTEGER(iwp) ::  l     !<
-    INTEGER(iwp) ::  m     !<
+    REAL(wp),DIMENSION(:,:) ::  a  !<
+    REAL(wp),DIMENSION(:,:) ::  b  !<
+    REAL(wp),DIMENSION(:,:) ::  c  !<
+    REAL(wp),DIMENSION(:,:) ::  d  !<
+
+    REAL(wp),DIMENSION(:),INTENT(IN) ::  trigs  !<
+
+    REAL(wp) ::  a0      !<
+    REAL(wp) ::  a1      !<
+    REAL(wp) ::  a10     !<
+    REAL(wp) ::  a11     !<
+    REAL(wp) ::  a2      !<
+    REAL(wp) ::  a20     !<
+    REAL(wp) ::  a21     !<
+    REAL(wp) ::  a3      !<
+    REAL(wp) ::  a4      !<
+    REAL(wp) ::  a5      !<
+    REAL(wp) ::  a6      !<
+    REAL(wp) ::  b0      !<
+    REAL(wp) ::  b1      !<
+    REAL(wp) ::  b10     !<
+    REAL(wp) ::  b11     !<
+    REAL(wp) ::  b2      !<
+    REAL(wp) ::  b20     !<
+    REAL(wp) ::  b21     !<
+    REAL(wp) ::  b3      !<
+    REAL(wp) ::  b4      !<
+    REAL(wp) ::  b5      !<
+    REAL(wp) ::  b6      !<
+    REAL(wp) ::  c1      !<
+    REAL(wp) ::  c2      !<
+    REAL(wp) ::  c3      !<
+    REAL(wp) ::  c4      !<
+    REAL(wp) ::  c5      !<
+    REAL(wp) ::  qrt5    !<
+    REAL(wp) ::  s1      !<
+    REAL(wp) ::  s2      !<
+    REAL(wp) ::  s3      !<
+    REAL(wp) ::  s4      !<
+    REAL(wp) ::  s5      !<
+    REAL(wp) ::  sin36   !<
+    REAL(wp) ::  sin45   !<
+    REAL(wp) ::  sin60   !<
+    REAL(wp) ::  sin72   !<
+    REAL(wp) ::  z       !<
+    REAL(wp) ::  zqrt5   !<
+    REAL(wp) ::  zsin36  !<
+    REAL(wp) ::  zsin45  !<
+    REAL(wp) ::  zsin60  !<
+    REAL(wp) ::  zsin72  !<
+
 
 
@@ -2531,6 +2583,6 @@
     iink  = la * inc1
     jink  = la * inc2
-    ijump = (ifac-1) * iink
-    kstop = ( n-ifac ) / ( 2*ifac )
+    ijump = ( ifac - 1 ) * iink
+    kstop = ( n - ifac ) / ( 2 * ifac )
 
     ibase = 0
@@ -2553,20 +2605,20 @@
           ib = ia + iink
           ja = 1
-          jb = ja + (2*m-la) * inc2
+          jb = ja + ( 2 * m - la ) * inc2
 
           IF ( la /= m )  THEN
 
              DO  l = 1, la
-                i = ibase
-                j = jbase
+                i         = ibase
+                j         = jbase
                 c(:,ja+j) = a(:,ia+i) + a(:,ib+i)
                 c(:,jb+j) = a(:,ia+i) - a(:,ib+i)
-                ibase = ibase + inc1
-                jbase = jbase + inc2
-             ENDDO
-
-             ja = ja + jink
-             jink = 2 * jink
-             jb = jb - jink
+                ibase     = ibase + inc1
+                jbase     = jbase + inc2
+             ENDDO
+
+             ja    = ja + jink
+             jink  = 2 * jink
+             jb    = jb - jink
              ibase = ibase + ijump
              ijump = 2 * ijump + iink
@@ -2582,8 +2634,8 @@
                       i = ibase
                       j = jbase
-                      c(:,ja+j) = a(:,ia+i) + (c1*a(:,ib+i)+s1*b(:,ib+i))
-                      c(:,jb+j) = a(:,ia+i) - (c1*a(:,ib+i)+s1*b(:,ib+i))
-                      d(:,ja+j) = (c1*b(:,ib+i)-s1*a(:,ib+i)) + b(:,ia+i)
-                      d(:,jb+j) = (c1*b(:,ib+i)-s1*a(:,ib+i)) - b(:,ia+i)
+                      c(:,ja+j) = a(:,ia+i) + ( c1 * a(:,ib+i) + s1 * b(:,ib+i) )
+                      c(:,jb+j) = a(:,ia+i) - ( c1 * a(:,ib+i) + s1 * b(:,ib+i) )
+                      d(:,ja+j) = ( c1 * b(:,ib+i) - s1 * a(:,ib+i) ) + b(:,ia+i)
+                      d(:,jb+j) = ( c1 * b(:,ib+i) - s1 * a(:,ib+i) ) - b(:,ia+i)
                       ibase = ibase + inc1
                       jbase = jbase + inc2
@@ -2611,10 +2663,10 @@
           ELSE
 
-             z = 1.0_wp/REAL(n,KIND=wp)
-             DO  l = 1, la
-                i = ibase
-                j = jbase
-                c(:,ja+j) = z*(a(:,ia+i)+a(:,ib+i))
-                c(:,jb+j) = z*(a(:,ia+i)-a(:,ib+i))
+             z = 1.0_wp / REAL( n, KIND = wp )
+             DO  l = 1, la
+                i = ibase
+                j = jbase
+                c(:,ja+j) = z * ( a(:,ia+i) + a(:,ib+i) )
+                c(:,jb+j) = z * ( a(:,ia+i) - a(:,ib+i) )
                 ibase = ibase + inc1
                 jbase = jbase + inc2
@@ -2631,5 +2683,5 @@
           ic = ib + iink
           ja = 1
-          jb = ja + (2*m-la) * inc2
+          jb = ja + ( 2 * m - la ) * inc2
           jc = jb
 
@@ -2637,17 +2689,17 @@
 
              DO  l = 1, la
-                i = ibase
-                j = jbase
-                c(:,ja+j) = a(:,ia+i) + (a(:,ib+i)+a(:,ic+i))
-                c(:,jb+j) = a(:,ia+i) - 0.5_wp*(a(:,ib+i)+a(:,ic+i))
-                d(:,jb+j) = sin60*(a(:,ic+i)-a(:,ib+i))
-                ibase = ibase + inc1
-                jbase = jbase + inc2
-             ENDDO
-
-             ja = ja + jink
-             jink = 2 * jink
-             jb = jb + jink
-             jc = jc - jink
+                i         = ibase
+                j         = jbase
+                c(:,ja+j) = a(:,ia+i) + ( a(:,ib+i) + a(:,ic+i) )
+                c(:,jb+j) = a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) )
+                d(:,jb+j) = sin60 * ( a(:,ic+i) - a(:,ib+i) )
+                ibase     = ibase + inc1
+                jbase     = jbase + inc2
+             ENDDO
+
+             ja    = ja + jink
+             jink  = 2 * jink
+             jb    = jb + jink
+             jc    = jc - jink
              ibase = ibase + ijump
              ijump = 2 * ijump + iink
@@ -2667,11 +2719,15 @@
                       i = ibase
                       j = jbase
-                      DO  ipl = 1, SIZE(a,1)
-                         a1 = (c1*a(ipl,ib+i)+s1*b(ipl,ib+i)) + (c2*a(ipl,ic+i)+s2*b(ipl,ic+i))
-                         b1 = (c1*b(ipl,ib+i)-s1*a(ipl,ib+i)) + (c2*b(ipl,ic+i)-s2*a(ipl,ic+i))
-                         a2 = a(ipl,ia+i) - 0.5_wp*a1
-                         b2 = b(ipl,ia+i) - 0.5_wp*b1
-                         a3 = sin60*((c1*a(ipl,ib+i)+s1*b(ipl,ib+i))-(c2*a(ipl,ic+i)+s2*b(ipl,ic+i)))
-                         b3 = sin60*((c1*b(ipl,ib+i)-s1*a(ipl,ib+i))-(c2*b(ipl,ic+i)-s2*a(ipl,ic+i)))
+                      DO  ipl = 1, SIZE( a, 1 )
+                         a1 = ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) )                              &
+                            + ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) )
+                         b1 = ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) )                              &
+                            + ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) )
+                         a2 = a(ipl,ia+i) - 0.5_wp * a1
+                         b2 = b(ipl,ia+i) - 0.5_wp * b1
+                         a3 = sin60 * ( ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) )                    &
+                                      - ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) ) )
+                         b3 = sin60 * ( ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) )                    &
+                                      - ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) ) )
                          c(ipl,ja+j) = a(ipl,ia+i) + a1
                          d(ipl,ja+j) = b(ipl,ia+i) + b1
@@ -2679,5 +2735,5 @@
                          d(ipl,jb+j) = b2 - a3
                          c(ipl,jc+j) = a2 - b3
-                         d(ipl,jc+j) = -(b2+a3)
+                         d(ipl,jc+j) = - ( b2 + a3 )
                       ENDDO
                       ibase = ibase + inc1
@@ -2685,7 +2741,7 @@
                    ENDDO
                    ibase = ibase + ijump
-                   ja = ja + jink
-                   jb = jb + jink
-                   jc = jc - jink
+                   ja    = ja + jink
+                   jb    = jb + jink
+                   jc    = jc - jink
                 ENDDO
 
@@ -2698,7 +2754,7 @@
                 i = ibase
                 j = jbase
-                c(:,ja+j) = a(:,ia+i) + 0.5_wp*(a(:,ib+i)-a(:,ic+i))
-                d(:,ja+j) = -sin60*(a(:,ib+i)+a(:,ic+i))
-                c(:,jb+j) = a(:,ia+i) - (a(:,ib+i)-a(:,ic+i))
+                c(:,ja+j) = a(:,ia+i) + 0.5_wp * ( a(:,ib+i) - a(:,ic+i) )
+                d(:,ja+j) = - sin60 * ( a(:,ib+i) + a(:,ic+i) )
+                c(:,jb+j) = a(:,ia+i) - ( a(:,ib+i) - a(:,ic+i) )
                 ibase = ibase + inc1
                 jbase = jbase + inc2
@@ -2707,12 +2763,12 @@
           ELSE
 
-             z = 1.0_wp / REAL( n, KIND=wp )
-             zsin60 = z*sin60
-             DO  l = 1, la
-                i = ibase
-                j = jbase
-                c(:,ja+j) = z*(a(:,ia+i)+(a(:,ib+i)+a(:,ic+i)))
-                c(:,jb+j) = z*(a(:,ia+i)-0.5_wp*(a(:,ib+i)+a(:,ic+i)))
-                d(:,jb+j) = zsin60*(a(:,ic+i)-a(:,ib+i))
+             z = 1.0_wp / REAL( n, KIND = wp )
+             zsin60 = z * sin60
+             DO  l = 1, la
+                i = ibase
+                j = jbase
+                c(:,ja+j) = z * ( a(:,ia+i) + ( a(:,ib+i) + a(:,ic+i) ) )
+                c(:,jb+j) = z * ( a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) )
+                d(:,jb+j) = zsin60 * ( a(:,ic+i) - a(:,ib+i) )
                 ibase = ibase + inc1
                 jbase = jbase + inc2
@@ -2730,6 +2786,6 @@
           id = ic + iink
           ja = 1
-          jb = ja + (2*m-la) * inc2
-          jc = jb + 2*m*inc2
+          jb = ja + ( 2 * m - la ) * inc2
+          jc = jb + 2 * m * inc2
           jd = jb
 
@@ -2737,9 +2793,8 @@
 
              DO  l = 1, la
-
-                i = ibase
-                j = jbase
-                c(:,ja+j) = (a(:,ia+i)+a(:,ic+i)) + (a(:,ib+i)+a(:,id+i))
-                c(:,jc+j) = (a(:,ia+i)+a(:,ic+i)) - (a(:,ib+i)+a(:,id+i))
+                i = ibase
+                j = jbase
+                c(:,ja+j) = ( a(:,ia+i) + a(:,ic+i) ) + ( a(:,ib+i) + a(:,id+i) )
+                c(:,jc+j) = ( a(:,ia+i) + a(:,ic+i) ) - ( a(:,ib+i) + a(:,id+i) )
                 c(:,jb+j) = a(:,ia+i) - a(:,ic+i)
                 d(:,jb+j) = a(:,id+i) - a(:,ib+i)
@@ -2748,9 +2803,9 @@
              ENDDO
 
-             ja = ja + jink
-             jink = 2 * jink
-             jb = jb + jink
-             jc = jc - jink
-             jd = jd - jink
+             ja    = ja + jink
+             jink  = 2 * jink
+             jb    = jb + jink
+             jc    = jc - jink
+             jd    = jd - jink
              ibase = ibase + ijump
              ijump = 2 * ijump + iink
@@ -2772,13 +2827,17 @@
                       i = ibase
                       j = jbase
-                      DO  ipl = 1, SIZE(a,1)
-                         a0 = a(ipl,ia+i) + (c2*a(ipl,ic+i)+s2*b(ipl,ic+i))
-                         a2 = a(ipl,ia+i) - (c2*a(ipl,ic+i)+s2*b(ipl,ic+i))
-                         a1 = (c1*a(ipl,ib+i)+s1*b(ipl,ib+i)) + (c3*a(ipl,id+i)+s3*b(ipl,id+i))
-                         a3 = (c1*a(ipl,ib+i)+s1*b(ipl,ib+i)) - (c3*a(ipl,id+i)+s3*b(ipl,id+i))
-                         b0 = b(ipl,ia+i) + (c2*b(ipl,ic+i)-s2*a(ipl,ic+i))
-                         b2 = b(ipl,ia+i) - (c2*b(ipl,ic+i)-s2*a(ipl,ic+i))
-                         b1 = (c1*b(ipl,ib+i)-s1*a(ipl,ib+i)) + (c3*b(ipl,id+i)-s3*a(ipl,id+i))
-                         b3 = (c1*b(ipl,ib+i)-s1*a(ipl,ib+i)) - (c3*b(ipl,id+i)-s3*a(ipl,id+i))
+                      DO  ipl = 1, SIZE( a, 1 )
+                         a0 = a(ipl,ia+i) + ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) )
+                         a2 = a(ipl,ia+i) - ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) )
+                         a1 = ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) )                              &
+                            + ( c3 * a(ipl,id+i) + s3 * b(ipl,id+i) )
+                         a3 = ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) )                              &
+                            - ( c3 * a(ipl,id+i) + s3 * b(ipl,id+i) )
+                         b0 = b(ipl,ia+i) + ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) )
+                         b2 = b(ipl,ia+i) - ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) )
+                         b1 = ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) )                              &
+                            + ( c3 * b(ipl,id+i) - s3 * a(ipl,id+i) )
+                         b3 = ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) )                              &
+                            - ( c3 * b(ipl,id+i) - s3 * a(ipl,id+i) )
                          c(ipl,ja+j) = a0 + a1
                          c(ipl,jc+j) = a0 - a1
@@ -2788,5 +2847,5 @@
                          c(ipl,jd+j) = a2 - b3
                          d(ipl,jb+j) = b2 - a3
-                         d(ipl,jd+j) = -(b2+a3)
+                         d(ipl,jd+j) = - ( b2 + a3 )
                       ENDDO
                       ibase = ibase + inc1
@@ -2809,8 +2868,8 @@
                 i = ibase
                 j = jbase
-                c(:,ja+j) = a(:,ia+i) + sin45*(a(:,ib+i)-a(:,id+i))
-                c(:,jb+j) = a(:,ia+i) - sin45*(a(:,ib+i)-a(:,id+i))
-                d(:,ja+j) = -a(:,ic+i) - sin45*(a(:,ib+i)+a(:,id+i))
-                d(:,jb+j) = a(:,ic+i) - sin45*(a(:,ib+i)+a(:,id+i))
+                c(:,ja+j) =   a(:,ia+i) + sin45 * ( a(:,ib+i) - a(:,id+i) )
+                c(:,jb+j) =   a(:,ia+i) - sin45 * ( a(:,ib+i) - a(:,id+i) )
+                d(:,ja+j) = - a(:,ic+i) - sin45 * ( a(:,ib+i) + a(:,id+i) )
+                d(:,jb+j) =   a(:,ic+i) - sin45 * ( a(:,ib+i) + a(:,id+i) )
                 ibase = ibase + inc1
                 jbase = jbase + inc2
@@ -2819,12 +2878,12 @@
           ELSE
 
-             z = 1.0_wp / REAL( n, KIND=wp )
-             DO  l = 1, la
-                i = ibase
-                j = jbase
-                c(:,ja+j) = z*((a(:,ia+i)+a(:,ic+i))+(a(:,ib+i)+a(:,id+i)))
-                c(:,jc+j) = z*((a(:,ia+i)+a(:,ic+i))-(a(:,ib+i)+a(:,id+i)))
-                c(:,jb+j) = z*(a(:,ia+i)-a(:,ic+i))
-                d(:,jb+j) = z*(a(:,id+i)-a(:,ib+i))
+             z = 1.0_wp / REAL( n, KIND = wp )
+             DO  l = 1, la
+                i = ibase
+                j = jbase
+                c(:,ja+j) = z * ( ( a(:,ia+i) + a(:,ic+i) ) + ( a(:,ib+i) + a(:,id+i) ) )
+                c(:,jc+j) = z * ( ( a(:,ia+i) + a(:,ic+i) ) - ( a(:,ib+i) + a(:,id+i) ) )
+                c(:,jb+j) = z * ( a(:,ia+i) - a(:,ic+i) )
+                d(:,jb+j) = z * ( a(:,id+i) - a(:,ib+i) )
                 ibase = ibase + inc1
                 jbase = jbase + inc2
@@ -2843,6 +2902,6 @@
           ie = id + iink
           ja = 1
-          jb = ja + (2*m-la) * inc2
-          jc = jb + 2*m*inc2
+          jb = ja + ( 2 * m - la ) * inc2
+          jc = jb + 2 * m * inc2
           jd = jc
           je = jb
@@ -2853,27 +2912,27 @@
                 i = ibase
                 j = jbase
-                DO  ipl = 1, SIZE(a,1)
+                DO  ipl = 1, SIZE( a, 1 )
                    a1 = a(ipl,ib+i) + a(ipl,ie+i)
                    a3 = a(ipl,ib+i) - a(ipl,ie+i)
                    a2 = a(ipl,ic+i) + a(ipl,id+i)
                    a4 = a(ipl,ic+i) - a(ipl,id+i)
-                   a5 = a(ipl,ia+i) - 0.25_wp*(a1+a2)
-                   a6 = qrt5*(a1-a2)
-                   c(ipl,ja+j) = a(ipl,ia+i) + (a1+a2)
+                   a5 = a(ipl,ia+i) - 0.25_wp * ( a1 + a2 )
+                   a6 = qrt5 * ( a1 - a2 )
+                   c(ipl,ja+j) = a(ipl,ia+i) + ( a1 + a2 )
                    c(ipl,jb+j) = a5 + a6
                    c(ipl,jc+j) = a5 - a6
-                   d(ipl,jb+j) = -sin72*a3 - sin36*a4
-                   d(ipl,jc+j) = -sin36*a3 + sin72*a4
-                ENDDO
-                ibase = ibase + inc1
-                jbase = jbase + inc2
-             ENDDO
-
-             ja = ja + jink
-             jink = 2 * jink
-             jb = jb + jink
-             jc = jc + jink
-             jd = jd - jink
-             je = je - jink
+                   d(ipl,jb+j) = - sin72 * a3 - sin36 * a4
+                   d(ipl,jc+j) = - sin36 * a3 + sin72 * a4
+                ENDDO
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+             ja    = ja + jink
+             jink  = 2 * jink
+             jb    = jb + jink
+             jc    = jc + jink
+             jd    = jd - jink
+             je    = je - jink
              ibase = ibase + ijump
              ijump = 2 * ijump + iink
@@ -2899,34 +2958,48 @@
                       j = jbase
                       DO  ipl = 1, SIZE(a,1)
-                         a1 = (c1*a(ipl,ib+i)+s1*b(ipl,ib+i)) + (c4*a(ipl,ie+i)+s4*b(ipl,ie+i))
-                         a3 = (c1*a(ipl,ib+i)+s1*b(ipl,ib+i)) - (c4*a(ipl,ie+i)+s4*b(ipl,ie+i))
-                         a2 = (c2*a(ipl,ic+i)+s2*b(ipl,ic+i)) + (c3*a(ipl,id+i)+s3*b(ipl,id+i))
-                         a4 = (c2*a(ipl,ic+i)+s2*b(ipl,ic+i)) - (c3*a(ipl,id+i)+s3*b(ipl,id+i))
-                         b1 = (c1*b(ipl,ib+i)-s1*a(ipl,ib+i)) + (c4*b(ipl,ie+i)-s4*a(ipl,ie+i))
-                         b3 = (c1*b(ipl,ib+i)-s1*a(ipl,ib+i)) - (c4*b(ipl,ie+i)-s4*a(ipl,ie+i))
-                         b2 = (c2*b(ipl,ic+i)-s2*a(ipl,ic+i)) + (c3*b(ipl,id+i)-s3*a(ipl,id+i))
-                         b4 = (c2*b(ipl,ic+i)-s2*a(ipl,ic+i)) - (c3*b(ipl,id+i)-s3*a(ipl,id+i))
-                         a5 = a(ipl,ia+i) - 0.25_wp*(a1+a2)
-                         a6 = qrt5*(a1-a2)
-                         b5 = b(ipl,ia+i) - 0.25_wp*(b1+b2)
-                         b6 = qrt5*(b1-b2)
+                         a1 = ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) ) + ( c4 * a(ipl,ie+i) + s4    &
+                              * b(ipl,ie+i) )
+                         a3 = ( c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i) ) - ( c4 * a(ipl,ie+i) + s4    &
+                              * b(ipl,ie+i) )
+                         a2 = ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) ) + ( c3 * a(ipl,id+i) + s3    &
+                              * b(ipl,id+i) )
+                         a4 = ( c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i) ) - ( c3 * a(ipl,id+i) + s3    &
+                              * b(ipl,id+i) )
+
+                         b1 = ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) ) + ( c4 * b(ipl,ie+i) - s4    &
+                              * a(ipl,ie+i) )
+                         b3 = ( c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i) ) - ( c4 * b(ipl,ie+i) - s4    &
+                              * a(ipl,ie+i) )
+                         b2 = ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) ) + ( c3 * b(ipl,id+i) - s3    &
+                              * a(ipl,id+i) )
+                         b4 = ( c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i) ) - ( c3 * b(ipl,id+i) - s3    &
+                              * a(ipl,id+i) )
+
+                         a5 = a(ipl,ia+i) - 0.25_wp * ( a1 + a2 )
+                         a6 = qrt5 * ( a1 - a2 )
+                         b5 = b(ipl,ia+i) - 0.25_wp * ( b1 + b2 )
+                         b6 = qrt5 * ( b1 - b2)
+
                          a10 = a5 + a6
                          a20 = a5 - a6
                          b10 = b5 + b6
                          b20 = b5 - b6
-                         a11 = sin72*b3 + sin36*b4
-                         a21 = sin36*b3 - sin72*b4
-                         b11 = sin72*a3 + sin36*a4
-                         b21 = sin36*a3 - sin72*a4
-                         c(ipl,ja+j) = a(ipl,ia+i) + (a1+a2)
+                         a11 = sin72 * b3 + sin36 * b4
+                         a21 = sin36 * b3 - sin72 * b4
+
+                         b11 = sin72 * a3 + sin36 * a4
+                         b21 = sin36 * a3 - sin72 * a4
+
+                         c(ipl,ja+j) = a(ipl,ia+i) + ( a1 + a2 )
                          c(ipl,jb+j) = a10 + a11
                          c(ipl,je+j) = a10 - a11
                          c(ipl,jc+j) = a20 + a21
                          c(ipl,jd+j) = a20 - a21
-                         d(ipl,ja+j) = b(ipl,ia+i) + (b1+b2)
+
+                         d(ipl,ja+j) = b(ipl,ia+i) + ( b1 + b2 )
                          d(ipl,jb+j) = b10 - b11
-                         d(ipl,je+j) = -(b10+b11)
+                         d(ipl,je+j) = - ( b10 + b11 )
                          d(ipl,jc+j) = b20 - b21
-                         d(ipl,jd+j) = -(b20+b21)
+                         d(ipl,jd+j) = - ( b20 + b21 )
                       ENDDO
                       ibase = ibase + inc1
@@ -2934,9 +3007,9 @@
                    ENDDO
                    ibase = ibase + ijump
-                   ja = ja + jink
-                   jb = jb + jink
-                   jc = jc + jink
-                   jd = jd - jink
-                   je = je - jink
+                   ja    = ja + jink
+                   jb    = jb + jink
+                   jc    = jc + jink
+                   jd    = jd - jink
+                   je    = je - jink
                 ENDDO
 
@@ -2949,16 +3022,18 @@
                 i = ibase
                 j = jbase
-                DO  ipl = 1, SIZE(a,1)
+                DO  ipl = 1, SIZE( a, 1 )
                    a1 = a(ipl,ib+i) + a(ipl,ie+i)
                    a3 = a(ipl,ib+i) - a(ipl,ie+i)
                    a2 = a(ipl,ic+i) + a(ipl,id+i)
                    a4 = a(ipl,ic+i) - a(ipl,id+i)
-                   a5 = a(ipl,ia+i) + 0.25_wp*(a3-a4)
-                   a6 = qrt5*(a3+a4)
+                   a5 = a(ipl,ia+i) + 0.25_wp * ( a3 - a4 )
+                   a6 = qrt5 * ( a3 + a4 )
+
                    c(ipl,ja+j) = a5 + a6
                    c(ipl,jb+j) = a5 - a6
-                   c(ipl,jc+j) = a(ipl,ia+i) - (a3-a4)
-                   d(ipl,ja+j) = -sin36*a1 - sin72*a2
-                   d(ipl,jb+j) = -sin72*a1 + sin36*a2
+                   c(ipl,jc+j) = a(ipl,ia+i) - ( a3 - a4 )
+
+                   d(ipl,ja+j) = - sin36 * a1 - sin72 * a2
+                   d(ipl,jb+j) = - sin72 * a1 + sin36 * a2
                 ENDDO
                 ibase = ibase + inc1
@@ -2968,5 +3043,5 @@
           ELSE
 
-             z = 1.0_wp / REAL( n, KIND=wp )
+             z = 1.0_wp / REAL( n, KIND = wp )
              zqrt5  = z * qrt5
              zsin36 = z * sin36
@@ -2975,16 +3050,18 @@
                 i = ibase
                 j = jbase
-                DO  ipl = 1, SIZE(a,1)
+                DO  ipl = 1, SIZE( a, 1 )
                    a1 = a(ipl,ib+i) + a(ipl,ie+i)
                    a3 = a(ipl,ib+i) - a(ipl,ie+i)
                    a2 = a(ipl,ic+i) + a(ipl,id+i)
                    a4 = a(ipl,ic+i) - a(ipl,id+i)
-                   a5 = z*(a(ipl,ia+i)-0.25_wp*(a1+a2))
-                   a6 = zqrt5*(a1-a2)
-                   c(ipl,ja+j) = z*(a(ipl,ia+i)+(a1+a2))
+                   a5 = z * ( a(ipl,ia+i) - 0.25_wp * ( a1 + a2 ) )
+                   a6 = zqrt5 * ( a1 - a2 )
+
+                   c(ipl,ja+j) = z * ( a(ipl,ia+i) + ( a1 + a2 ) )
                    c(ipl,jb+j) = a5 + a6
                    c(ipl,jc+j) = a5 - a6
-                   d(ipl,jb+j) = -zsin72*a3 - zsin36*a4
-                   d(ipl,jc+j) = -zsin36*a3 + zsin72*a4
+
+                   d(ipl,jb+j) = - zsin72 * a3 - zsin36 * a4
+                   d(ipl,jc+j) = - zsin36 * a3 + zsin72 * a4
                 ENDDO
                 ibase = ibase + inc1
@@ -3005,7 +3082,7 @@
           if = ie + iink
           ja = 1
-          jb = ja + (2*m-la) * inc2
-          jc = jb + 2*m*inc2
-          jd = jc + 2*m*inc2
+          jb = ja + ( 2 * m - la ) * inc2
+          jc = jb + 2 * m * inc2
+          jd = jc + 2 * m * inc2
           je = jc
           jf = jb
@@ -3016,13 +3093,15 @@
                 i = ibase
                 j = jbase
-                DO  ipl = 1, SIZE(a,1)
-                   a11 = (a(ipl,ic+i)+a(ipl,if+i)) + (a(ipl,ib+i)+a(ipl,ie+i))
-                   c(ipl,ja+j) = (a(ipl,ia+i)+a(ipl,id+i)) + a11
-                   c(ipl,jc+j) = (a(ipl,ia+i)+a(ipl,id+i)-0.5_wp*a11)
-                   d(ipl,jc+j) = sin60*((a(ipl,ic+i)+a(ipl,if+i))-(a(ipl,ib+i)+a(ipl,ie+i)))
-                   a11 = (a(ipl,ic+i)-a(ipl,if+i)) + (a(ipl,ie+i)-a(ipl,ib+i))
-                   c(ipl,jb+j) = (a(ipl,ia+i)-a(ipl,id+i)) - 0.5_wp*a11
-                   d(ipl,jb+j) = sin60*((a(ipl,ie+i)-a(ipl,ib+i))-(a(ipl,ic+i)-a(ipl,if+i)))
-                   c(ipl,jd+j) = (a(ipl,ia+i)-a(ipl,id+i)) + a11
+                DO  ipl = 1, SIZE( a, 1 )
+                   a11         = ( a(ipl,ic+i) + a(ipl,if+i) ) + ( a(ipl,ib+i) + a(ipl,ie+i) )
+                   c(ipl,ja+j) = ( a(ipl,ia+i) + a(ipl,id+i) ) + a11
+                   c(ipl,jc+j) = ( a(ipl,ia+i) + a(ipl,id+i) - 0.5_wp * a11 )
+                   d(ipl,jc+j) = sin60 * ( ( a(ipl,ic+i) + a(ipl,if+i) ) -                         &
+                                 ( a(ipl,ib+i) + a(ipl,ie+i) ) )
+                   a11         = ( a(ipl,ic+i) - a(ipl,if+i) ) + ( a(ipl,ie+i) - a(ipl,ib+i) )
+                   c(ipl,jb+j) = ( a(ipl,ia+i) - a(ipl,id+i) ) - 0.5_wp * a11
+                   d(ipl,jb+j) = sin60 * ( ( a(ipl,ie+i) - a(ipl,ib+i) ) -                         &
+                                 ( a(ipl,ic+i) - a(ipl,if+i) ) )
+                   c(ipl,jd+j) = ( a(ipl,ia+i) - a(ipl,id+i) ) + a11
                 END DO
                 ibase = ibase + inc1
@@ -3061,37 +3140,37 @@
                       i = ibase
                       j = jbase
-                      DO  ipl = 1, SIZE(a,1)
-                         a1 = c1*a(ipl,ib+i) + s1*b(ipl,ib+i)
-                         b1 = c1*b(ipl,ib+i) - s1*a(ipl,ib+i)
-                         a2 = c2*a(ipl,ic+i) + s2*b(ipl,ic+i)
-                         b2 = c2*b(ipl,ic+i) - s2*a(ipl,ic+i)
-                         a3 = c3*a(ipl,id+i) + s3*b(ipl,id+i)
-                         b3 = c3*b(ipl,id+i) - s3*a(ipl,id+i)
-                         a4 = c4*a(ipl,ie+i) + s4*b(ipl,ie+i)
-                         b4 = c4*b(ipl,ie+i) - s4*a(ipl,ie+i)
-                         a5 = c5*a(ipl,if+i) + s5*b(ipl,if+i)
-                         b5 = c5*b(ipl,if+i) - s5*a(ipl,if+i)
-                         a11 = (a2+a5) + (a1+a4)
-                         a20 = (a(ipl,ia+i)+a3) - 0.5_wp*a11
-                         a21 = sin60*((a2+a5)-(a1+a4))
-                         b11 = (b2+b5) + (b1+b4)
-                         b20 = (b(ipl,ia+i)+b3) - 0.5_wp*b11
-                         b21 = sin60*((b2+b5)-(b1+b4))
-                         c(ipl,ja+j) = (a(ipl,ia+i)+a3) + a11
-                         d(ipl,ja+j) = (b(ipl,ia+i)+b3) + b11
+                      DO  ipl = 1, SIZE( a, 1 )
+                         a1 = c1 * a(ipl,ib+i) + s1 * b(ipl,ib+i)
+                         b1 = c1 * b(ipl,ib+i) - s1 * a(ipl,ib+i)
+                         a2 = c2 * a(ipl,ic+i) + s2 * b(ipl,ic+i)
+                         b2 = c2 * b(ipl,ic+i) - s2 * a(ipl,ic+i)
+                         a3 = c3 * a(ipl,id+i) + s3 * b(ipl,id+i)
+                         b3 = c3 * b(ipl,id+i) - s3 * a(ipl,id+i)
+                         a4 = c4 * a(ipl,ie+i) + s4 * b(ipl,ie+i)
+                         b4 = c4 * b(ipl,ie+i) - s4 * a(ipl,ie+i)
+                         a5 = c5 * a(ipl,if+i) + s5 * b(ipl,if+i)
+                         b5 = c5 * b(ipl,if+i) - s5 * a(ipl,if+i)
+                         a11 = ( a2 + a5 ) + ( a1 + a4 )
+                         a20 = ( a(ipl,ia+i) + a3 ) - 0.5_wp * a11
+                         a21 = sin60 * ( ( a2 + a5 ) - ( a1 + a4 ) )
+                         b11 = ( b2 + b5 ) + ( b1 + b4 )
+                         b20 = ( b(ipl,ia+i) + b3 ) - 0.5_wp * b11
+                         b21 = sin60 * ( ( b2 + b5 ) - ( b1 + b4 ) )
+                         c(ipl,ja+j) = ( a(ipl,ia+i) + a3 ) + a11
+                         d(ipl,ja+j) = ( b(ipl,ia+i) + b3 ) + b11
                          c(ipl,jc+j) = a20 - b21
                          d(ipl,jc+j) = a21 + b20
                          c(ipl,je+j) = a20 + b21
                          d(ipl,je+j) = a21 - b20
-                         a11 = (a2-a5) + (a4-a1)
-                         a20 = (a(ipl,ia+i)-a3) - 0.5_wp*a11
-                         a21 = sin60*((a4-a1)-(a2-a5))
-                         b11 = (b5-b2) - (b4-b1)
-                         b20 = (b3-b(ipl,ia+i)) - 0.5_wp*b11
-                         b21 = sin60*((b5-b2)+(b4-b1))
+                         a11 = ( a2 - a5 ) + ( a4 - a1 )
+                         a20 = ( a(ipl,ia+i) - a3 ) - 0.5_wp * a11
+                         a21 = sin60 * ( ( a4 - a1 ) - ( a2 - a5 ) )
+                         b11 = ( b5 - b2 ) - ( b4 - b1 )
+                         b20 = ( b3 - b(ipl,ia+i) ) - 0.5_wp * b11
+                         b21 = sin60 * ( ( b5 - b2 ) + ( b4 - b1 ) )
                          c(ipl,jb+j) = a20 - b21
                          d(ipl,jb+j) = a21 - b20
-                         c(ipl,jd+j) = a11 + (a(ipl,ia+i)-a3)
-                         d(ipl,jd+j) = b11 + (b3-b(ipl,ia+i))
+                         c(ipl,jd+j) = a11 + ( a(ipl,ia+i) - a3 )
+                         d(ipl,jd+j) = b11 + ( b3 - b(ipl,ia+i) )
                          c(ipl,jf+j) = a20 + b21
                          d(ipl,jf+j) = a21 + b20
@@ -3117,10 +3196,14 @@
                 i = ibase
                 j = jbase
-                c(:,ja+j) = (a(:,ia+i)+0.5_wp*(a(:,ic+i)-a(:,ie+i))) + sin60*(a(:,ib+i)-a(:,if+i))
-                d(:,ja+j) = -(a(:,id+i)+0.5_wp*(a(:,ib+i)+a(:,if+i))) - sin60*(a(:,ic+i)+a(:,ie+i))
-                c(:,jb+j) = a(:,ia+i) - (a(:,ic+i)-a(:,ie+i))
-                d(:,jb+j) = a(:,id+i) - (a(:,ib+i)+a(:,if+i))
-                c(:,jc+j) = (a(:,ia+i)+0.5_wp*(a(:,ic+i)-a(:,ie+i))) - sin60*(a(:,ib+i)-a(:,if+i))
-                d(:,jc+j) = -(a(:,id+i)+0.5_wp*(a(:,ib+i)+a(:,if+i))) + sin60*(a(:,ic+i)+a(:,ie+i))
+                c(:,ja+j) = ( a(:,ia+i) + 0.5_wp * ( a(:,ic+i) - a(:,ie+i) ) ) + sin60             &
+                          * ( a(:,ib+i) - a(:,if+i) )
+                d(:,ja+j) = - ( a(:,id+i) + 0.5_wp * ( a(:,ib+i) + a(:,if+i) ) ) - sin60           &
+                            * ( a(:,ic+i) + a(:,ie+i) )
+                c(:,jb+j) = a(:,ia+i) - ( a(:,ic+i) - a(:,ie+i) )
+                d(:,jb+j) = a(:,id+i) - ( a(:,ib+i) + a(:,if+i) )
+                c(:,jc+j) = ( a(:,ia+i) + 0.5_wp * ( a(:,ic+i) - a(:,ie+i) ) ) - sin60             &
+                          * ( a(:,ib+i) - a(:,if+i) )
+                d(:,jc+j) = - ( a(:,id+i) + 0.5_wp * ( a(:,ib+i) + a(:,if+i) ) ) + sin60           &
+                            * ( a(:,ic+i) + a(:,ie+i) )
                 ibase = ibase + inc1
                 jbase = jbase + inc2
@@ -3129,18 +3212,20 @@
           ELSE
 
-             z = 1.0_wp/REAL(n,KIND=wp)
-             zsin60 = z*sin60
-             DO  l = 1, la
-                i = ibase
-                j = jbase
-                DO  ipl = 1, SIZE(a,1)
-                   a11 = (a(ipl,ic+i)+a(ipl,if+i)) + (a(ipl,ib+i)+a(ipl,ie+i))
-                   c(ipl,ja+j) = z*((a(ipl,ia+i)+a(ipl,id+i))+a11)
-                   c(ipl,jc+j) = z*((a(ipl,ia+i)+a(ipl,id+i))-0.5_wp*a11)
-                   d(ipl,jc+j) = zsin60*((a(ipl,ic+i)+a(ipl,if+i))-(a(ipl,ib+i)+a(ipl,ie+i)))
-                   a11 = (a(ipl,ic+i)-a(ipl,if+i)) + (a(ipl,ie+i)-a(ipl,ib+i))
-                   c(ipl,jb+j) = z*((a(ipl,ia+i)-a(ipl,id+i))-0.5_wp*a11)
-                   d(ipl,jb+j) = zsin60*((a(ipl,ie+i)-a(ipl,ib+i))-(a(ipl,ic+i)-a(ipl,if+i)))
-                   c(ipl,jd+j) = z*((a(ipl,ia+i)-a(ipl,id+i))+a11)
+             z = 1.0_wp / REAL( n, KIND = wp )
+             zsin60 = z * sin60
+             DO  l = 1, la
+                i = ibase
+                j = jbase
+                DO  ipl = 1, SIZE( a, 1 )
+                   a11 = ( a (ipl,ic+i) + a(ipl,if+i) ) + ( a(ipl,ib+i) + a(ipl,ie+i) )
+                   c(ipl,ja+j) = z * ( ( a(ipl,ia+i) + a(ipl,id+i) ) + a11 )
+                   c(ipl,jc+j) = z * ( ( a(ipl,ia+i) + a(ipl,id+i) ) - 0.5_wp * a11 )
+                   d(ipl,jc+j) = zsin60 * ( ( a(ipl,ic+i) + a(ipl,if+i) )                          &
+                                          - ( a(ipl,ib+i) + a(ipl,ie+i) ) )
+                   a11 = ( a(ipl,ic+i) - a(ipl,if+i) )  + ( a(ipl,ie+i) - a(ipl,ib+i) )
+                   c(ipl,jb+j) = z * ( ( a(ipl,ia+i) - a(ipl,id+i) ) - 0.5_wp * a11 )
+                   d(ipl,jb+j) = zsin60 * ( ( a(ipl,ie+i) - a(ipl,ib+i) )                          &
+                                          - ( a(ipl,ic+i) - a(ipl,if+i) ) )
+                   c(ipl,jd+j) = z * ( ( a(ipl,ia+i) - a(ipl,id+i) ) + a11 )
                 ENDDO
                 ibase = ibase + inc1
@@ -3169,8 +3254,8 @@
           ja = 1
           jb = ja + la * inc2
-          jc = jb + 2*m*inc2
-          jd = jc + 2*m*inc2
-          je = jd + 2*m*inc2
-          z = 1.0_wp / REAL( n, KIND=wp )
+          jc = jb + 2 * m * inc2
+          jd = jc + 2 * m * inc2
+          je = jd + 2 * m * inc2
+          z = 1.0_wp / REAL( n, KIND = wp )
           zsin45 = z * SQRT( 0.5_wp )
 
@@ -3178,12 +3263,18 @@
              i = ibase
              j = jbase
-             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))))
-             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))))
-             c(:,jc+j) = z*((a(:,ia+i)+a(:,ie+i))-(a(:,ic+i)+a(:,ig+i)))
-             d(:,jc+j) = z*((a(:,id+i)+a(:,ih+i))-(a(:,ib+i)+a(:,if+i)))
-             c(:,jb+j) = z*(a(:,ia+i)-a(:,ie+i)) + zsin45*((a(:,ih+i)-a(:,id+i))-(a(:,if+i)-a(:,ib+i)))
-             c(:,jd+j) = z*(a(:,ia+i)-a(:,ie+i)) - zsin45*((a(:,ih+i)-a(:,id+i))-(a(:,if+i)-a(:,ib+i)))
-             d(:,jb+j) = zsin45*((a(:,ih+i)-a(:,id+i))+(a(:,if+i)-a(:,ib+i))) + z*(a(:,ig+i)-a(:,ic+i))
-             d(:,jd+j) = zsin45*((a(:,ih+i)-a(:,id+i))+(a(:,if+i)-a(:,ib+i))) - z*(a(:,ig+i)-a(:,ic+i))
+             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) ) ) )
+             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) ) ) )
+             c(:,jc+j) = z * ( ( a(:,ia+i) + a(:,ie+i) ) - ( a(:,ic+i) + a(:,ig+i) ) )
+             d(:,jc+j) = z * ( ( a(:,id+i) + a(:,ih+i) ) - ( a(:,ib+i) + a(:,if+i) ) )
+             c(:,jb+j) = z * ( a(:,ia+i) - a(:,ie+i) ) + zsin45 * ( ( a(:,ih+i) - a(:,id+i) )      &
+                         - ( a(:,if+i) - a(:,ib+i) ) )
+             c(:,jd+j) = z * ( a(:,ia+i) - a(:,ie+i) ) - zsin45 * ( ( a(:,ih+i) - a(:,id+i) )      &
+                         - ( a(:,if+i) - a(:,ib+i) ) )
+             d(:,jb+j) = zsin45 * ( ( a(:,ih+i) - a(:,id+i) ) + ( a(:,if+i) - a(:,ib+i) ) )        &
+                         + z * ( a(:,ig+i) - a(:,ic+i) )
+             d(:,jd+j) = zsin45 * ( ( a(:,ih+i) - a(:,id+i) ) + ( a(:,if+i) - a(:,ib+i) ) )        &
+                         - z * ( a(:,ig+i) - a(:,ic+i) )
              ibase = ibase + inc1
              jbase = jbase + inc2
@@ -3194,9 +3285,9 @@
  END SUBROUTINE qpassm_vec
 
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
 ! Description:
 ! ------------
 !> Same as qpassm, but for backward fft
-!------------------------------------------------------------------------------!
+!--------------------------------------------------------------------------------------------------!
  SUBROUTINE rpassm_vec(a, b, c, d, trigs, inc1, inc2, n, ifac, la, ierr )
 
@@ -3205,80 +3296,80 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) ::  ierr !<
-    INTEGER(iwp) ::  ifac !<
-    INTEGER(iwp) ::  inc1 !<
-    INTEGER(iwp) ::  inc2 !<
-    INTEGER(iwp) ::  la   !<
-    INTEGER(iwp) ::  n    !<
+
+    INTEGER(iwp) ::  i      !<
+    INTEGER(iwp) ::  ia     !<
+    INTEGER(iwp) ::  ib     !<
+    INTEGER(iwp) ::  ibase  !<
+    INTEGER(iwp) ::  ic     !<
+    INTEGER(iwp) ::  id     !<
+    INTEGER(iwp) ::  ie     !<
+    INTEGER(iwp) ::  ierr   !<
+    INTEGER(iwp) ::  ifac   !<
+    INTEGER(iwp) ::  if     !<
+    INTEGER(iwp) ::  igo    !<
+    INTEGER(iwp) ::  iink   !<
+    INTEGER(iwp) ::  inc1   !<
+    INTEGER(iwp) ::  inc2   !<
+    INTEGER(iwp) ::  ipl    !<  loop index parallel loop
+    INTEGER(iwp) ::  j      !<
+    INTEGER(iwp) ::  ja     !<
+    INTEGER(iwp) ::  jb     !<
+    INTEGER(iwp) ::  jbase  !<
+    INTEGER(iwp) ::  jc     !<
+    INTEGER(iwp) ::  jd     !<
+    INTEGER(iwp) ::  je     !<
+    INTEGER(iwp) ::  jf     !<
+    INTEGER(iwp) ::  jg     !<
+    INTEGER(iwp) ::  jh     !<
+    INTEGER(iwp) ::  jink   !<
+    INTEGER(iwp) ::  jump   !<
+    INTEGER(iwp) ::  k      !<
+    INTEGER(iwp) ::  kb     !<
+    INTEGER(iwp) ::  kc     !<
+    INTEGER(iwp) ::  kd     !<
+    INTEGER(iwp) ::  ke     !<
+    INTEGER(iwp) ::  kf     !<
+    INTEGER(iwp) ::  kstop  !<
+    INTEGER(iwp) ::  l      !<
+    INTEGER(iwp) ::  la     !<
+    INTEGER(iwp) ::  m      !<
+    INTEGER(iwp) ::  n      !<
 !
 !-- Arrays are dimensioned with n
-    REAL(wp),DIMENSION(:,:) ::  a     !<
-    REAL(wp),DIMENSION(:,:) ::  b     !<
-    REAL(wp),DIMENSION(:,:) ::  c     !<
-    REAL(wp),DIMENSION(:,:) ::  d     !<
-    REAL(wp),DIMENSION(:),INTENT(IN) ::  trigs !<
-
-    REAL(wp) ::  c1     !<
-    REAL(wp) ::  c2     !<
-    REAL(wp) ::  c3     !<
-    REAL(wp) ::  c4     !<
-    REAL(wp) ::  c5     !<
-    REAL(wp) ::  qqrt5  !<
-    REAL(wp) ::  qrt5   !<
-    REAL(wp) ::  s1     !<
-    REAL(wp) ::  s2     !<
-    REAL(wp) ::  s3     !<
-    REAL(wp) ::  s4     !<
-    REAL(wp) ::  s5     !<
-    REAL(wp) ::  sin36  !<
-    REAL(wp) ::  sin45  !<
-    REAL(wp) ::  sin60  !<
-    REAL(wp) ::  sin72  !<
-    REAL(wp) ::  ssin36 !<
-    REAL(wp) ::  ssin45 !<
-    REAL(wp) ::  ssin60 !<
-    REAL(wp) ::  ssin72 !<
-
-    INTEGER(iwp) ::  i     !<
-    INTEGER(iwp) ::  ia    !<
-    INTEGER(iwp) ::  ib    !<
-    INTEGER(iwp) ::  ibase !<
-    INTEGER(iwp) ::  ic    !<
-    INTEGER(iwp) ::  id    !<
-    INTEGER(iwp) ::  ie    !<
-    INTEGER(iwp) ::  if    !<
-    INTEGER(iwp) ::  igo   !<
-    INTEGER(iwp) ::  iink  !<
-    INTEGER(iwp) ::  ipl   !<  loop index parallel loop
-    INTEGER(iwp) ::  j     !<
-    INTEGER(iwp) ::  ja    !<
-    INTEGER(iwp) ::  jb    !<
-    INTEGER(iwp) ::  jbase !<
-    INTEGER(iwp) ::  jc    !<
-    INTEGER(iwp) ::  jd    !<
-    INTEGER(iwp) ::  je    !<
-    INTEGER(iwp) ::  jf    !<
-    INTEGER(iwp) ::  jg    !<
-    INTEGER(iwp) ::  jh    !<
-    INTEGER(iwp) ::  jink  !<
-    INTEGER(iwp) ::  jump  !<
-    INTEGER(iwp) ::  k     !<
-    INTEGER(iwp) ::  kb    !<
-    INTEGER(iwp) ::  kc    !<
-    INTEGER(iwp) ::  kd    !<
-    INTEGER(iwp) ::  ke    !<
-    INTEGER(iwp) ::  kf    !<
-    INTEGER(iwp) ::  kstop !<
-    INTEGER(iwp) ::  l     !<
-    INTEGER(iwp) ::  m     !<
-
-    REAL(wp) ::  a10       !<
-    REAL(wp) ::  a11       !<
-    REAL(wp) ::  a20       !<
-    REAL(wp) ::  a21       !<
-    REAL(wp) ::  b10       !<
-    REAL(wp) ::  b11       !<
-    REAL(wp) ::  b20       !<
-    REAL(wp) ::  b21       !<
+    REAL(wp),DIMENSION(:,:) ::  a  !<
+    REAL(wp),DIMENSION(:,:) ::  b  !<
+    REAL(wp),DIMENSION(:,:) ::  c  !<
+    REAL(wp),DIMENSION(:,:) ::  d  !<
+    REAL(wp),DIMENSION(:),INTENT(IN) ::  trigs  !<
+
+    REAL(wp) ::  a10     !<
+    REAL(wp) ::  a11     !<
+    REAL(wp) ::  a20     !<
+    REAL(wp) ::  a21     !<
+    REAL(wp) ::  b10     !<
+    REAL(wp) ::  b11     !<
+    REAL(wp) ::  b20     !<
+    REAL(wp) ::  b21     !<
+    REAL(wp) ::  c1      !<
+    REAL(wp) ::  c2      !<
+    REAL(wp) ::  c3      !<
+    REAL(wp) ::  c4      !<
+    REAL(wp) ::  c5      !<
+    REAL(wp) ::  qqrt5   !<
+    REAL(wp) ::  qrt5    !<
+    REAL(wp) ::  s1      !<
+    REAL(wp) ::  s2      !<
+    REAL(wp) ::  s3      !<
+    REAL(wp) ::  s4      !<
+    REAL(wp) ::  s5      !<
+    REAL(wp) ::  sin36   !<
+    REAL(wp) ::  sin45   !<
+    REAL(wp) ::  sin60   !<
+    REAL(wp) ::  sin72   !<
+    REAL(wp) ::  ssin36  !<
+    REAL(wp) ::  ssin45  !<
+    REAL(wp) ::  ssin60  !<
+    REAL(wp) ::  ssin72  !<
+
 
 
@@ -3292,6 +3383,6 @@
     iink = la * inc1
     jink = la * inc2
-    jump = (ifac-1) * jink
-    kstop = (n-ifac) / (2*ifac)
+    jump = ( ifac - 1 ) * jink
+    kstop = ( n - ifac ) / ( 2 * ifac )
 
     ibase = 0
@@ -3311,5 +3402,5 @@
 
           ia = 1
-          ib = ia + (2*m-la) * inc1
+          ib = ia + ( 2 * m - la ) * inc1
           ja = 1
           jb = ja + jink
@@ -3345,6 +3436,6 @@
                       c(:,ja+j) = a(:,ia+i) + a(:,ib+i)
                       d(:,ja+j) = b(:,ia+i) - b(:,ib+i)
-                      c(:,jb+j) = c1*(a(:,ia+i)-a(:,ib+i)) - s1*(b(:,ia+i)+b(:,ib+i))
-                      d(:,jb+j) = s1*(a(:,ia+i)-a(:,ib+i)) + c1*(b(:,ia+i)+b(:,ib+i))
+                      c(:,jb+j) = c1 * ( a(:,ia+i) - a(:,ib+i) ) - s1 * ( b(:,ia+i) + b(:,ib+i) )
+                      d(:,jb+j) = s1 * ( a(:,ia+i) - a(:,ib+i) ) + c1 * ( b(:,ia+i) + b(:,ib+i) )
                       ibase = ibase + inc1
                       jbase = jbase + inc2
@@ -3364,5 +3455,5 @@
                 j = jbase
                 c(:,ja+j) = a(:,ia+i)
-                c(:,jb+j) = -b(:,ia+i)
+                c(:,jb+j) = - b(:,ia+i)
                 ibase = ibase + inc1
                 jbase = jbase + inc2
@@ -3374,6 +3465,6 @@
                 i = ibase
                 j = jbase
-                c(:,ja+j) = 2.0_wp*(a(:,ia+i)+a(:,ib+i))
-                c(:,jb+j) = 2.0_wp*(a(:,ia+i)-a(:,ib+i))
+                c(:,ja+j) = 2.0_wp * ( a(:,ia+i) + a(:,ib+i) )
+                c(:,jb+j) = 2.0_wp * ( a(:,ia+i) - a(:,ib+i) )
                 ibase = ibase + inc1
                 jbase = jbase + inc2
@@ -3387,5 +3478,5 @@
 
           ia = 1
-          ib = ia + (2*m-la) * inc1
+          ib = ia + ( 2 * m - la ) * inc1
           ic = ib
           ja = 1
@@ -3399,6 +3490,6 @@
                 j = jbase
                 c(:,ja+j) = a(:,ia+i) + a(:,ib+i)
-                c(:,jb+j) = (a(:,ia+i)-0.5_wp*a(:,ib+i)) - (sin60*(b(:,ib+i)))
-                c(:,jc+j) = (a(:,ia+i)-0.5_wp*a(:,ib+i)) + (sin60*(b(:,ib+i)))
+                c(:,jb+j) = ( a(:,ia+i) - 0.5_wp * a(:,ib+i) ) - ( sin60 * ( b(:,ib+i) ) )
+                c(:,jc+j) = ( a(:,ia+i) - 0.5_wp * a(:,ib+i) ) + ( sin60 * ( b(:,ib+i) ) )
                 ibase = ibase + inc1
                 jbase = jbase + inc2
@@ -3424,14 +3515,22 @@
                       i = ibase
                       j = jbase
-                      c(:,ja+j) = a(:,ia+i) + (a(:,ib+i)+a(:,ic+i))
-                      d(:,ja+j) = b(:,ia+i) + (b(:,ib+i)-b(:,ic+i))
-                      c(:,jb+j) = c1*((a(:,ia+i)-0.5_wp*(a(:,ib+i)+a(:,ic+i)))-(sin60*(b(:,ib+i)+ b(:,ic+i)))) &
-                                - s1*((b(:,ia+i)-0.5_wp*(b(:,ib+i)-b(:,ic+i)))+(sin60*(a(:,ib+i)- a(:,ic+i))))
-                      d(:,jb+j) = s1*((a(:,ia+i)-0.5_wp*(a(:,ib+i)+a(:,ic+i)))-(sin60*(b(:,ib+i)+ b(:,ic+i)))) &
-                                + c1*((b(:,ia+i)-0.5_wp*(b(:,ib+i)-b(:,ic+i)))+(sin60*(a(:,ib+i)- a(:,ic+i))))
-                      c(:,jc+j) = c2*((a(:,ia+i)-0.5_wp*(a(:,ib+i)+a(:,ic+i)))+(sin60*(b(:,ib+i)+ b(:,ic+i)))) &
-                                - s2*((b(:,ia+i)-0.5_wp*(b(:,ib+i)-b(:,ic+i)))-(sin60*(a(:,ib+i)- a(:,ic+i))))
-                      d(:,jc+j) = s2*((a(:,ia+i)-0.5_wp*(a(:,ib+i)+a(:,ic+i)))+(sin60*(b(:,ib+i)+ b(:,ic+i)))) &
-                                + c2*((b(:,ia+i)-0.5_wp*(b(:,ib+i)-b(:,ic+i)))-(sin60*(a(:,ib+i)- a(:,ic+i))))
+                      c(:,ja+j) = a(:,ia+i) + ( a(:,ib+i) + a(:,ic+i) )
+                      d(:,ja+j) = b(:,ia+i) + ( b(:,ib+i) - b(:,ic+i) )
+                      c(:,jb+j) = c1 * ( ( a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) )        &
+                                  - ( sin60 * ( b(:,ib+i) + b(:,ic+i) ) ) ) - s1 * ( ( b(:,ia+i)   &
+                                  - 0.5_wp * ( b(:,ib+i) - b(:,ic+i) ) ) + ( sin60 * ( a(:,ib+i)   &
+                                  - a(:,ic+i) ) ) )
+                      d(:,jb+j) = s1 * ( ( a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) )        &
+                                  - ( sin60 * ( b(:,ib+i) + b(:,ic+i) ) ) ) + c1 * ( ( b(:,ia+i)   &
+                                  - 0.5_wp * ( b(:,ib+i) - b(:,ic+i) ) ) + ( sin60 * ( a(:,ib+i)   &
+                                  - a(:,ic+i) ) ) )
+                      c(:,jc+j) = c2 * ( ( a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) )        &
+                                  + ( sin60 * ( b(:,ib+i) + b(:,ic+i) ) ) ) - s2 * ( ( b(:,ia+i)   &
+                                  - 0.5_wp * ( b(:,ib+i) - b(:,ic+i) ) ) - ( sin60 * ( a(:,ib+i)   &
+                                  - a(:,ic+i) ) ) )
+                      d(:,jc+j) = s2 * ( ( a(:,ia+i) - 0.5_wp * ( a(:,ib+i) + a(:,ic+i) ) )        &
+                                  + ( sin60 * ( b(:,ib+i) + b(:,ic+i) ) ) ) + c2 * ( ( b(:,ia+i)   &
+                                  - 0.5_wp * ( b(:,ib+i) - b(:,ic+i) ) ) - ( sin60 * ( a(:,ib+i)   &
+                                  - a(:,ic+i) ) ) )
                       ibase = ibase + inc1
                       jbase = jbase + inc2
@@ -3452,6 +3551,6 @@
                 j = jbase
                 c(:,ja+j) = a(:,ia+i) + a(:,ib+i)
-                c(:,jb+j) = (0.5_wp*a(:,ia+i)-a(:,ib+i)) - (sin60*b(:,ia+i))
-                c(:,jc+j) = -(0.5_wp*a(:,ia+i)-a(:,ib+i)) - (sin60*b(:,ia+i))
+                c(:,jb+j) = ( 0.5_wp * a(:,ia+i) - a(:,ib+i) ) - ( sin60 * b(:,ia+i) )
+                c(:,jc+j) = - ( 0.5_wp * a(:,ia+i) - a(:,ib+i) ) - ( sin60 * b(:,ia+i) )
                 ibase = ibase + inc1
                 jbase = jbase + inc2
@@ -3464,7 +3563,7 @@
                 i = ibase
                 j = jbase
-                c(:,ja+j) = 2.0_wp*(a(:,ia+i)+a(:,ib+i))
-                c(:,jb+j) = (2.0_wp*a(:,ia+i)-a(:,ib+i)) - (ssin60*b(:,ib+i))
-                c(:,jc+j) = (2.0_wp*a(:,ia+i)-a(:,ib+i)) + (ssin60*b(:,ib+i))
+                c(:,ja+j) = 2.0_wp * ( a(:,ia+i) + a(:,ib+i) )
+                c(:,jb+j) = ( 2.0_wp * a(:,ia+i) - a(:,ib+i) ) - ( ssin60 * b(:,ib+i) )
+                c(:,jc+j) = ( 2.0_wp * a(:,ia+i) - a(:,ib+i) ) + ( ssin60 * b(:,ib+i) )
                 ibase = ibase + inc1
                 jbase = jbase + inc2
@@ -3478,6 +3577,6 @@
 
           ia = 1
-          ib = ia + (2*m-la) * inc1
-          ic = ib + 2*m*inc1
+          ib = ia + ( 2 * m - la ) * inc1
+          ic = ib + 2 * m * inc1
           id = ib
           ja = 1
@@ -3491,8 +3590,8 @@
                 i = ibase
                 j = jbase
-                c(:,ja+j) = (a(:,ia+i)+a(:,ic+i)) + a(:,ib+i)
-                c(:,jb+j) = (a(:,ia+i)-a(:,ic+i)) - b(:,ib+i)
-                c(:,jc+j) = (a(:,ia+i)+a(:,ic+i)) - a(:,ib+i)
-                c(:,jd+j) = (a(:,ia+i)-a(:,ic+i)) + b(:,ib+i)
+                c(:,ja+j) = ( a(:,ia+i) + a(:,ic+i) ) + a(:,ib+i)
+                c(:,jb+j) = ( a(:,ia+i) - a(:,ic+i) ) - b(:,ib+i)
+                c(:,jc+j) = ( a(:,ia+i) + a(:,ic+i) ) - a(:,ib+i)
+                c(:,jd+j) = ( a(:,ia+i) - a(:,ic+i) ) + b(:,ib+i)
                 ibase = ibase + inc1
                 jbase = jbase + inc2
@@ -3522,18 +3621,18 @@
                       i = ibase
                       j = jbase
-                      c(:,ja+j) = (a(:,ia+i)+a(:,ic+i)) + (a(:,ib+i)+a(:,id+i))
-                      d(:,ja+j) = (b(:,ia+i)-b(:,ic+i)) + (b(:,ib+i)-b(:,id+i))
-                      c(:,jc+j) = c2*((a(:,ia+i)+a(:,ic+i))-(a(:,ib+i)+a(:,id+i))) - s2*((b(:,ia+i)-b(:,ic+i))&
-                                -(b(:,ib+i)-b(:,id+i)))
-                      d(:,jc+j) = s2*((a(:,ia+i)+a(:,ic+i))-(a(:,ib+i)+a(:,id+i))) + c2*((b(:,ia+i)-b(:,ic+i))&
-                                -(b(:,ib+i)-b(:,id+i)))
-                      c(:,jb+j) = c1*((a(:,ia+i)-a(:,ic+i))-(b(:,ib+i)+b(:,id+i))) - s1*((b(:,ia+i)+b(:,ic+i))&
-                                +(a(:,ib+i)-a(:,id+i)))
-                      d(:,jb+j) = s1*((a(:,ia+i)-a(:,ic+i))-(b(:,ib+i)+b(:,id+i))) + c1*((b(:,ia+i)+b(:,ic+i))&
-                                +(a(:,ib+i)-a(:,id+i)))
-                      c(:,jd+j) = c3*((a(:,ia+i)-a(:,ic+i))+(b(:,ib+i)+b(:,id+i))) - s3*((b(:,ia+i)+b(:,ic+i))&
-                                -(a(:,ib+i)-a(:,id+i)))
-                      d(:,jd+j) = s3*((a(:,ia+i)-a(:,ic+i))+(b(:,ib+i)+b(:,id+i))) + c3*((b(:,ia+i)+b(:,ic+i))&
-                                -(a(:,ib+i)-a(:,id+i)))
+                      c(:,ja+j) = ( a(:,ia+i) + a(:,ic+i) ) + ( a(:,ib+i) + a(:,id+i) )
+                      d(:,ja+j) = ( b(:,ia+i) - b(:,ic+i) ) + ( b(:,ib+i) - b(:,id+i) )
+                      c(:,jc+j) = c2 * ( ( a(:,ia+i) + a(:,ic+i) ) - ( a(:,ib+i) + a(:,id+i) ) )   &
+                                  - s2 * ( ( b(:,ia+i) - b(:,ic+i) ) - ( b(:,ib+i) - b(:,id+i) ) )
+                      d(:,jc+j) = s2 * ( ( a(:,ia+i) + a(:,ic+i) ) - ( a(:,ib+i) + a(:,id+i) ) )   &
+                                  + c2 * ( ( b(:,ia+i) - b(:,ic+i) ) - ( b(:,ib+i) - b(:,id+i) ) )
+                      c(:,jb+j) = c1 * ( ( a(:,ia+i) - a(:,ic+i) ) - ( b(:,ib+i) + b(:,id+i) ) )   &
+                                  - s1 * ( ( b(:,ia+i) + b(:,ic+i) ) + ( a(:,ib+i) - a(:,id+i) ) )
+                      d(:,jb+j) = s1 * ( ( a(:,ia+i) - a(:,ic+i) ) - ( b(:,ib+i) + b(:,id+i) ) )   &
+                                  + c1 * ( ( b(:,ia+i) + b(:,ic+i) ) + ( a(:,ib+i) - a(:,id+i) ) )
+                      c(:,jd+j) = c3 * ( ( a(:,ia+i) - a(:,ic+i) ) + ( b(:,ib+i) + b(:,id+i) ) )   &
+                                  - s3 * ( ( b(:,ia+i) + b(:,ic+i) ) - ( a(:,ib+i) - a(:,id+i) ) )
+                      d(:,jd+j) = s3 * ( ( a(:,ia+i) - a(:,ic+i) ) + ( b(:,ib+i) + b(:,id+i) ) )   &
+                                  + c3 * ( ( b(:,ia+i) + b(:,ic+i) ) - ( a(:,ib+i) - a(:,id+i) ) )
                       ibase = ibase + inc1
                       jbase = jbase + inc2
@@ -3556,7 +3655,7 @@
                 j = jbase
                 c(:,ja+j) = a(:,ia+i) + a(:,ib+i)
-                c(:,jb+j) = sin45*((a(:,ia+i)-a(:,ib+i))-(b(:,ia+i)+b(:,ib+i)))
+                c(:,jb+j) = sin45 * ( ( a(:,ia+i) - a(:,ib+i) ) - ( b(:,ia+i) + b(:,ib+i) ) )
                 c(:,jc+j) = b(:,ib+i) - b(:,ia+i)
-                c(:,jd+j) = -sin45*((a(:,ia+i)-a(:,ib+i))+(b(:,ia+i)+b(:,ib+i)))
+                c(:,jd+j) = - sin45 * ( ( a(:,ia+i) - a(:,ib+i) ) + ( b(:,ia+i) + b(:,ib+i) ) )
                 ibase = ibase + inc1
                 jbase = jbase + inc2
@@ -3568,8 +3667,8 @@
                 i = ibase
                 j = jbase
-                c(:,ja+j) = 2.0_wp*((a(:,ia+i)+a(:,ic+i))+a(:,ib+i))
-                c(:,jb+j) = 2.0_wp*((a(:,ia+i)-a(:,ic+i))-b(:,ib+i))
-                c(:,jc+j) = 2.0_wp*((a(:,ia+i)+a(:,ic+i))-a(:,ib+i))
-                c(:,jd+j) = 2.0_wp*((a(:,ia+i)-a(:,ic+i))+b(:,ib+i))
+                c(:,ja+j) = 2.0_wp * ( ( a(:,ia+i) + a(:,ic+i) ) + a(:,ib+i) )
+                c(:,jb+j) = 2.0_wp * ( ( a(:,ia+i) - a(:,ic+i) ) - b(:,ib+i) )
+                c(:,jc+j) = 2.0_wp * ( ( a(:,ia+i) + a(:,ic+i) ) - a(:,ib+i) )
+                c(:,jd+j) = 2.0_wp * ( ( a(:,ia+i) - a(:,ic+i) ) + b(:,ib+i) )
                 ibase = ibase + inc1
                 jbase = jbase + inc2
@@ -3583,6 +3682,6 @@
 
           ia = 1
-          ib = ia + (2*m-la) * inc1
-          ic = ib + 2*m*inc1
+          ib = ia + ( 2 * m - la ) * inc1
+          ic = ib + 2 * m * inc1
           id = ic
           ie = ib
@@ -3598,13 +3697,13 @@
                 i = ibase
                 j = jbase
-                c(:,ja+j) = a(:,ia+i) + (a(:,ib+i)+a(:,ic+i))
-                c(:,jb+j) = ((a(:,ia+i)-0.25_wp*(a(:,ib+i)+a(:,ic+i)))+qrt5*(a(:,ib+i)-a(:,ic+i))) -        &
-                          (sin72*b(:,ib+i)+sin36*b(:,ic+i))
-                c(:,jc+j) = ((a(:,ia+i)-0.25_wp*(a(:,ib+i)+a(:,ic+i)))-qrt5*(a(:,ib+i)-a(:,ic+i))) -        &
-                          (sin36*b(:,ib+i)-sin72*b(:,ic+i))
-                c(:,jd+j) = ((a(:,ia+i)-0.25_wp*(a(:,ib+i)+a(:,ic+i)))-qrt5*(a(:,ib+i)-a(:,ic+i))) +        &
-                          (sin36*b(:,ib+i)-sin72*b(:,ic+i))
-                c(:,je+j) = ((a(:,ia+i)-0.25_wp*(a(:,ib+i)+a(:,ic+i)))+qrt5*(a(:,ib+i)-a(:,ic+i))) +        &
-                          (sin72*b(:,ib+i)+sin36*b(:,ic+i))
+                c(:,ja+j) = a(:,ia+i) + ( a(:,ib+i) + a(:,ic+i) )
+                c(:,jb+j) = ( ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) ) + qrt5           &
+                          * ( a(:,ib+i) - a(:,ic+i) ) ) - ( sin72 * b(:,ib+i) + sin36 * b(:,ic+i) )
+                c(:,jc+j) = ( ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) ) - qrt5           &
+                          * ( a(:,ib+i) - a(:,ic+i) ) ) - ( sin36 * b(:,ib+i) - sin72 * b(:,ic+i) )
+                c(:,jd+j) = ( ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) ) - qrt5           &
+                          * ( a(:,ib+i) - a(:,ic+i) ) ) + ( sin36 * b(:,ib+i) - sin72 * b(:,ic+i) )
+                c(:,je+j) = ( ( a(:,ia+i) - 0.25_wp * ( a(:,ib+i) + a(:,ic+i) ) ) + qrt5           &
+                          * ( a(:,ib+i) - a(:,ic+i) ) ) + ( sin72 * b(:,ib+i) + sin36 * b(:,ic+i) )
                 ibase = ibase + inc1
                 jbase = jbase + inc2
@@ -3639,27 +3738,37 @@
                       j = jbase
 !DIR$ IVDEP
-                      DO  ipl = 1, SIZE(a,1)
-                         a10      = (a(ipl,ia+i)-0.25_wp*((a(ipl,ib+i)+a(ipl,ie+i))+(a(ipl,ic+i)+a(ipl,id+i)))) +      &
-                                    qrt5*((a(ipl,ib+i)+a(ipl,ie+i))-(a(ipl,ic+i)+a(ipl,id+i)))
-                         a20      = (a(ipl,ia+i)-0.25_wp*((a(ipl,ib+i)+a(ipl,ie+i))+(a(ipl,ic+i)+a(ipl,id+i)))) -      &
-                                    qrt5*((a(ipl,ib+i)+a(ipl,ie+i))-(a(ipl,ic+i)+a(ipl,id+i)))
-                         b10      = (b(ipl,ia+i)-0.25_wp*((b(ipl,ib+i)-b(ipl,ie+i))+(b(ipl,ic+i)-b(ipl,id+i)))) +      &
-                                    qrt5*((b(ipl,ib+i)-b(ipl,ie+i))-(b(ipl,ic+i)-b(ipl,id+i)))
-                         b20      = (b(ipl,ia+i)-0.25_wp*((b(ipl,ib+i)-b(ipl,ie+i))+(b(ipl,ic+i)-b(ipl,id+i)))) -      &
-                                    qrt5*((b(ipl,ib+i)-b(ipl,ie+i))-(b(ipl,ic+i)-b(ipl,id+i)))
-                         a11      = sin72*(b(ipl,ib+i)+b(ipl,ie+i)) + sin36*(b(ipl,ic+i)+b(ipl,id+i))
-                         a21      = sin36*(b(ipl,ib+i)+b(ipl,ie+i)) - sin72*(b(ipl,ic+i)+b(ipl,id+i))
-                         b11      = sin72*(a(ipl,ib+i)-a(ipl,ie+i)) + sin36*(a(ipl,ic+i)-a(ipl,id+i))
-                         b21      = sin36*(a(ipl,ib+i)-a(ipl,ie+i)) - sin72*(a(ipl,ic+i)-a(ipl,id+i))
-                         c(ipl,ja+j) = a(ipl,ia+i) + ((a(ipl,ib+i)+a(ipl,ie+i))+(a(ipl,ic+i)+a(ipl,id+i)))
-                         d(ipl,ja+j) = b(ipl,ia+i) + ((b(ipl,ib+i)-b(ipl,ie+i))+(b(ipl,ic+i)-b(ipl,id+i)))
-                         c(ipl,jb+j) = c1*(a10     -a11     ) - s1*(b10     +b11     )
-                         d(ipl,jb+j) = s1*(a10     -a11     ) + c1*(b10     +b11     )
-                         c(ipl,je+j) = c4*(a10     +a11     ) - s4*(b10     -b11     )
-                         d(ipl,je+j) = s4*(a10     +a11     ) + c4*(b10     -b11     )
-                         c(ipl,jc+j) = c2*(a20     -a21     ) - s2*(b20     +b21     )
-                         d(ipl,jc+j) = s2*(a20     -a21     ) + c2*(b20     +b21     )
-                         c(ipl,jd+j) = c3*(a20     +a21     ) - s3*(b20     -b21     )
-                         d(ipl,jd+j) = s3*(a20     +a21     ) + c3*(b20     -b21     )
+                      DO  ipl = 1, SIZE( a, 1 )
+                         a10 = ( a (ipl,ia+i) - 0.25_wp * ( ( a(ipl,ib+i) + a(ipl,ie+i) )          &
+                             + ( a(ipl,ic+i) + a(ipl,id+i) ) ) ) + qrt5 * ( ( a(ipl,ib+i)          &
+                             + a(ipl,ie+i) ) - ( a(ipl,ic+i) + a(ipl,id+i) ) )
+                         a20 = ( a(ipl,ia+i) - 0.25_wp * ( ( a(ipl,ib+i) + a(ipl,ie+i) )           &
+                             + ( a(ipl,ic+i) + a(ipl,id+i) ) ) ) - qrt5 * ( ( a(ipl,ib+i)          &
+                             + a(ipl,ie+i) ) - ( a(ipl,ic+i) + a(ipl,id+i) ) )
+                         b10 = ( b(ipl,ia+i) - 0.25_wp * ( ( b(ipl,ib+i) - b(ipl,ie+i) )           &
+                             + ( b(ipl,ic+i) - b(ipl,id+i) ) ) ) + qrt5 * ( ( b(ipl,ib+i)          &
+                             - b(ipl,ie+i) ) - ( b(ipl,ic+i) - b(ipl,id+i) ) )
+                         b20 = ( b(ipl,ia+i) - 0.25_wp * ( ( b(ipl,ib+i) - b(ipl,ie+i) )           &
+                             + ( b(ipl,ic+i) - b(ipl,id+i) ) ) ) - qrt5 * ( ( b(ipl,ib+i)          &
+                             - b(ipl,ie+i) ) - ( b(ipl,ic+i) - b(ipl,id+i) ) )
+                         a11 = sin72 * ( b(ipl,ib+i) + b(ipl,ie+i) ) + sin36                       &
+                               * ( b(ipl,ic+i) + b(ipl,id+i) )
+                         a21 = sin36 * ( b(ipl,ib+i) + b(ipl,ie+i) ) - sin72                       &
+                               * ( b(ipl,ic+i) + b(ipl,id+i) )
+                         b11 = sin72 * ( a(ipl,ib+i) - a(ipl,ie+i) ) + sin36                       &
+                               * ( a(ipl,ic+i) - a(ipl,id+i) )
+                         b21 = sin36 * ( a(ipl,ib+i) - a(ipl,ie+i) ) - sin72                       &
+                               * ( a(ipl,ic+i) - a(ipl,id+i) )
+                         c(ipl,ja+j) = a(ipl,ia+i) + ( ( a(ipl,ib+i) + a(ipl,ie+i) )               &
+                                       + ( a(ipl,ic+i) + a(ipl,id+i) ) )
+                         d(ipl,ja+j) = b(ipl,ia+i) + ( ( b(ipl,ib+i) - b(ipl,ie+i) )               &
+                                       + ( b(ipl,ic+i) - b(ipl,id+i) ) )
+                         c(ipl,jb+j) = c1 * ( a10 - a11 ) - s1 * ( b10 + b11 )
+                         d(ipl,jb+j) = s1 * ( a10 - a11 ) + c1 * ( b10 + b11 )
+                         c(ipl,je+j) = c4 * ( a10 + a11 ) - s4 * ( b10 - b11 )
+                         d(ipl,je+j) = s4 * ( a10 + a11 ) + c4 * ( b10 - b11 )
+                         c(ipl,jc+j) = c2 * ( a20 - a21 ) - s2 * ( b20 + b21 )
+                         d(ipl,jc+j) = s2 * ( a20 - a21 ) + c2 * ( b20 + b21 )
+                         c(ipl,jd+j) = c3 * ( a20 + a21 ) - s3 * ( b20 - b21 )
+                         d(ipl,jd+j) = s3 * ( a20 + a21 ) + c3 * ( b20 - b21 )
                       ENDDO
                       ibase = ibase + inc1
@@ -3682,13 +3791,13 @@
                 i = ibase
                 j = jbase
-                c(:,ja+j) = (a(:,ia+i)+a(:,ib+i)) + a(:,ic+i)
-                c(:,jb+j) = (qrt5*(a(:,ia+i)-a(:,ib+i))+(0.25_wp*(a(:,ia+i)+a(:,ib+i))-a(:,ic+i))) -        &
-                          (sin36*b(:,ia+i)+sin72*b(:,ib+i))
-                c(:,je+j) = -(qrt5*(a(:,ia+i)-a(:,ib+i))+(0.25_wp*(a(:,ia+i)+a(:,ib+i))-a(:,ic+i))) -       &
-                          (sin36*b(:,ia+i)+sin72*b(:,ib+i))
-                c(:,jc+j) = (qrt5*(a(:,ia+i)-a(:,ib+i))-(0.25_wp*(a(:,ia+i)+a(:,ib+i))-a(:,ic+i))) -        &
-                          (sin72*b(:,ia+i)-sin36*b(:,ib+i))
-                c(:,jd+j) = -(qrt5*(a(:,ia+i)-a(:,ib+i))-(0.25_wp*(a(:,ia+i)+a(:,ib+i))-a(:,ic+i))) -       &
-                          (sin72*b(:,ia+i)-sin36*b(:,ib+i))
+                c(:,ja+j) = ( a(:,ia+i) + a(:,ib+i) ) + a(:,ic+i)
+                c(:,jb+j) = ( qrt5 * ( a(:,ia+i) - a(:,ib+i) ) + ( 0.25_wp * ( a(:,ia+i)           &
+                            + a(:,ib+i) ) - a(:,ic+i) ) ) - ( sin36 * b(:,ia+i) + sin72 * b(:,ib+i) )
+                c(:,je+j) = - ( qrt5 * ( a(:,ia+i) - a(:,ib+i) ) + ( 0.25_wp * ( a(:,ia+i)         &
+                            + a(:,ib+i) ) - a(:,ic+i) ) ) - ( sin36 * b(:,ia+i) + sin72 * b(:,ib+i) )
+                c(:,jc+j) = ( qrt5 * ( a(:,ia+i) - a(:,ib+i) ) - ( 0.25_wp * ( a(:,ia+i)           &
+                            + a(:,ib+i) ) - a(:,ic+i) ) ) - ( sin72 * b(:,ia+i) - sin36 * b(:,ib+i) )
+                c(:,jd+j) = - ( qrt5 * ( a(:,ia+i) - a(:,ib+i) ) - ( 0.25_wp * ( a(:,ia+i)         &
+                            + a(:,ib+i) ) - a(:,ic+i) ) ) - ( sin72 * b(:,ia+i) - sin36 * b(:,ib+i) )
                 ibase = ibase + inc1
                 jbase = jbase + inc2
@@ -3703,13 +3812,17 @@
                 i = ibase
                 j = jbase
-                c(:,ja+j) = 2.0_wp*(a(:,ia+i)+(a(:,ib+i)+a(:,ic+i)))
-                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)))  &
-                          - (ssin72*b(:,ib+i)+ssin36*b(:,ic+i))
-                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)))  &
-                          - (ssin36*b(:,ib+i)-ssin72*b(:,ic+i))
-                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)))  &
-                          + (ssin36*b(:,ib+i)-ssin72*b(:,ic+i))
-                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)))  &
-                          + (ssin72*b(:,ib+i)+ssin36*b(:,ic+i))
+                c(:,ja+j) = 2.0_wp * ( a(:,ia+i) + ( a(:,ib+i) + a(:,ic+i) ) )
+                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) ) )  - ( ssin72 * b(:,ib+i)          &
+                            + ssin36 * b(:,ic+i) )
+                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) ) ) - ( ssin36 * b(:,ib+i)           &
+                            - ssin72 * b(:,ic+i) )
+                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) ) ) + ( ssin36 * b(:,ib+i)           &
+                            - ssin72 * b(:,ic+i) )
+                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) ) ) + ( ssin72 * b(:,ib+i)           &
+                            + ssin36 * b(:,ic+i) )
                 ibase = ibase + inc1
                 jbase = jbase + inc2
@@ -3723,7 +3836,7 @@
 
           ia = 1
-          ib = ia + (2*m-la) * inc1
-          ic = ib + 2*m*inc1
-          id = ic + 2*m*inc1
+          ib = ia + ( 2 * m - la ) * inc1
+          ic = ib + 2 * m * inc1
+          id = ic + 2 * m * inc1
           ie = ic
           if = ib
@@ -3740,10 +3853,14 @@
                 i = ibase
                 j = jbase
-                c(:,ja+j) = (a(:,ia+i)+a(:,id+i)) + (a(:,ib+i)+a(:,ic+i))
-                c(:,jd+j) = (a(:,ia+i)-a(:,id+i)) - (a(:,ib+i)-a(:,ic+i))
-                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)))
-                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)))
-                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)))
-                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)))
+                c(:,ja+j) = ( a(:,ia+i) + a(:,id+i) ) +  ( a(:,ib+i) + a(:,ic+i) )
+                c(:,jd+j) = ( a(:,ia+i) - a(:,id+i) ) -  ( a(:,ib+i) - a(:,ic+i) )
+                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) ) )
+                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) ) )
+                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) ) )
+                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) ) )
                 ibase = ibase + inc1
                 jbase = jbase + inc2
@@ -3781,32 +3898,38 @@
                       i = ibase
                       j = jbase
-                      DO  ipl = 1, SIZE(a,1)
-                         a11      = (a(ipl,ie+i)+a(ipl,ib+i)) + (a(ipl,ic+i)+a(ipl,if+i))
-                         a20      = (a(ipl,ia+i)+a(ipl,id+i)) - 0.5_wp*a11
-                         a21      = sin60*((a(ipl,ie+i)+a(ipl,ib+i))-(a(ipl,ic+i)+a(ipl,if+i)))
-                         b11      = (b(ipl,ib+i)-b(ipl,ie+i)) + (b(ipl,ic+i)-b(ipl,if+i))
-                         b20      = (b(ipl,ia+i)-b(ipl,id+i)) - 0.5_wp*b11
-                         b21      = sin60*((b(ipl,ib+i)-b(ipl,ie+i))-(b(ipl,ic+i)-b(ipl,if+i)))
-
-                         c(ipl,ja+j) = (a(ipl,ia+i)+a(ipl,id+i)) + a11
-                         d(ipl,ja+j) = (b(ipl,ia+i)-b(ipl,id+i)) + b11
-                         c(ipl,jc+j) = c2*(a20     -b21     ) - s2*(b20     +a21     )
-                         d(ipl,jc+j) = s2*(a20     -b21     ) + c2*(b20     +a21     )
-                         c(ipl,je+j) = c4*(a20     +b21     ) - s4*(b20     -a21     )
-                         d(ipl,je+j) = s4*(a20     +b21     ) + c4*(b20     -a21     )
-
-                         a11      = (a(ipl,ie+i)-a(ipl,ib+i)) + (a(ipl,ic+i)-a(ipl,if+i))
-                         b11      = (b(ipl,ie+i)+b(ipl,ib+i)) - (b(ipl,ic+i)+b(ipl,if+i))
-                         a20      = (a(ipl,ia+i)-a(ipl,id+i)) - 0.5_wp*a11
-                         a21      = sin60*((a(ipl,ie+i)-a(ipl,ib+i))-(a(ipl,ic+i)-a(ipl,if+i)))
-                         b20      = (b(ipl,ia+i)+b(ipl,id+i)) + 0.5_wp*b11
-                         b21      = sin60*((b(ipl,ie+i)+b(ipl,ib+i))+(b(ipl,ic+i)+b(ipl,if+i)))
-
-                         c(ipl,jd+j) = c3*((a(ipl,ia+i)-a(ipl,id+i))+a11     ) - s3*((b(ipl,ia+i)+b(ipl,id+i))-b11     )
-                         d(ipl,jd+j) = s3*((a(ipl,ia+i)-a(ipl,id+i))+a11     ) + c3*((b(ipl,ia+i)+b(ipl,id+i))-b11     )
-                         c(ipl,jb+j) = c1*(a20     -b21     ) - s1*(b20     -a21     )
-                         d(ipl,jb+j) = s1*(a20     -b21     ) + c1*(b20     -a21     )
-                         c(ipl,jf+j) = c5*(a20     +b21     ) - s5*(b20     +a21     )
-                         d(ipl,jf+j) = s5*(a20     +b21     ) + c5*(b20     +a21     )
+                      DO  ipl = 1, SIZE( a, 1)
+                         a11 = ( a(ipl,ie+i) + a(ipl,ib+i) ) + ( a(ipl,ic+i) + a(ipl,if+i) )
+                         a20 = ( a(ipl,ia+i) + a(ipl,id+i) ) - 0.5_wp * a11
+                         a21 = sin60 * ( ( a(ipl,ie+i) + a(ipl,ib+i) )                             &
+                               - ( a(ipl,ic+i) + a(ipl,if+i) ) )
+                         b11 = ( b(ipl,ib+i) - b(ipl,ie+i) ) + ( b(ipl,ic+i) - b(ipl,if+i) )
+                         b20 = ( b(ipl,ia+i) - b(ipl,id+i) ) - 0.5_wp * b11
+                         b21 = sin60 * ( ( b(ipl,ib+i) - b(ipl,ie+i) )                             &
+                               - ( b(ipl,ic+i) - b(ipl,if+i) ) )
+
+                         c(ipl,ja+j) = ( a(ipl,ia+i) + a(ipl,id+i) ) + a11
+                         d(ipl,ja+j) = ( b(ipl,ia+i) - b(ipl,id+i) ) + b11
+                         c(ipl,jc+j) = c2 * ( a20 - b21 ) - s2 * ( b20 + a21 )
+                         d(ipl,jc+j) = s2 * ( a20 - b21 ) + c2 * ( b20 + a21 )
+                         c(ipl,je+j) = c4 * ( a20 + b21 ) - s4 * ( b20 - a21 )
+                         d(ipl,je+j) = s4 * ( a20 + b21 ) + c4 * ( b20 - a21 )
+
+                         a11 = ( a(ipl,ie+i) - a(ipl,ib+i) ) + ( a(ipl,ic+i) - a(ipl,if+i) )
+                         b11 = ( b(ipl,ie+i) + b(ipl,ib+i) ) - ( b(ipl,ic+i) + b(ipl,if+i) )
+                         a20 = ( a(ipl,ia+i) - a(ipl,id+i) ) - 0.5_wp * a11
+                         a21 = sin60 * ( ( a(ipl,ie+i) - a(ipl,ib+i) ) - ( a(ipl,ic+i)             &
+                               - a(ipl,if+i) ) )
+                         b20 = ( b(ipl,ia+i) + b(ipl,id+i) ) + 0.5_wp * b11
+                         b21 = sin60 * ( ( b(ipl,ie+i) + b(ipl,ib+i) ) + ( b(ipl,ic+i)             &
+                               + b(ipl,if+i) ) )
+
+                         c(ipl,jd+j) = c3 * ( ( a(ipl,ia+i) - a(ipl,id+i) ) + a11 )                &
+                                     - s3 * ( ( b(ipl,ia+i) + b(ipl,id+i) ) - b11 )
+                         d(ipl,jd+j) = s3 * ( ( a(ipl,ia+i) - a(ipl,id+i) ) + a11 )                &
+                                     + c3 * ( ( b(ipl,ia+i) + b(ipl,id+i) ) - b11 )
+                         c(ipl,jb+j) = c1 * ( a20 - b21 ) - s1 * ( b20 - a21 )
+                         d(ipl,jb+j) = s1 * ( a20 - b21 ) + c1 * ( b20 - a21 )
+                         c(ipl,jf+j) = c5 * ( a20 + b21 ) - s5 * ( b20 + a21 )
+                         d(ipl,jf+j) = s5 * ( a20 + b21 ) + c5 * ( b20 + a21 )
 
                       ENDDO
@@ -3831,11 +3954,15 @@
                 i = ibase
                 j = jbase
-                DO  ipl = 1, SIZE(a,1)
-                   c(ipl,ja+j) = a(ipl,ib+i) + (a(ipl,ia+i)+a(ipl,ic+i))
-                   c(ipl,jd+j) = b(ipl,ib+i) - (b(ipl,ia+i)+b(ipl,ic+i))
-                   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))
-                   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))
-                   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))
-                   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))
+                DO  ipl = 1, SIZE( a ,1 )
+                   c(ipl,ja+j) = a(ipl,ib+i) + ( a(ipl,ia+i) + a(ipl,ic+i) )
+                   c(ipl,jd+j) = b(ipl,ib+i) - ( b(ipl,ia+i) + b(ipl,ic+i) )
+                   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) )
+                   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) )
+                   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) )
+                   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) )
                 ENDDO
                 ibase = ibase + inc1
@@ -3849,10 +3976,16 @@
                 i = ibase
                 j = jbase
-                c(:,ja+j) = (2.0_wp*(a(:,ia+i)+a(:,id+i))) + (2.0_wp*(a(:,ib+i)+a(:,ic+i)))
-                c(:,jd+j) = (2.0_wp*(a(:,ia+i)-a(:,id+i))) - (2.0_wp*(a(:,ib+i)-a(:,ic+i)))
-                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)))
-                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)))
-                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)))
-                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)))
+                c(:,ja+j) = ( 2.0_wp * ( a(:,ia+i) + a(:,id+i) ) )                                 &
+                          + ( 2.0_wp * ( a(:,ib+i) + a(:,ic+i) ) )
+                c(:,jd+j) = ( 2.0_wp * ( a(:,ia+i) - a(:,id+i) ) )                                 &
+                          - ( 2.0_wp * ( a(:,ib+i) - a(:,ic+i) ) )
+                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) ) )
+                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) ) )
+                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) ) )
+                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) ) )
                 ibase = ibase + inc1
                 jbase = jbase + inc2
@@ -3871,8 +4004,8 @@
 
           ia = 1
-          ib = ia + la*inc1
-          ic = ib + 2*la*inc1
-          id = ic + 2*la*inc1
-          ie = id + 2*la*inc1
+          ib = ia + la * inc1
+          ic = ib + 2 * la * inc1
+          id = ic + 2 * la * inc1
+          ie = id + 2 * la * inc1
           ja = 1
           jb = ja + jink
@@ -3888,12 +4021,20 @@
              i = ibase
              j = jbase
-             c(:,ja+j) = 2.0_wp*(((a(:,ia+i)+a(:,ie+i))+a(:,ic+i))+(a(:,ib+i)+a(:,id+i)))
-             c(:,je+j) = 2.0_wp*(((a(:,ia+i)+a(:,ie+i))+a(:,ic+i))-(a(:,ib+i)+a(:,id+i)))
-             c(:,jc+j) = 2.0_wp*(((a(:,ia+i)+a(:,ie+i))-a(:,ic+i))-(b(:,ib+i)-b(:,id+i)))
-             c(:,jg+j) = 2.0_wp*(((a(:,ia+i)+a(:,ie+i))-a(:,ic+i))+(b(:,ib+i)-b(:,id+i)))
-             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)))
-             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)))
-             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)))
-             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)))
+             c(:,ja+j) = 2.0_wp * ( ( ( a(:,ia+i) + a(:,ie+i) ) + a(:,ic+i) )                      &
+                                                  + ( a(:,ib+i) + a(:,id+i) ) )
+             c(:,je+j) = 2.0_wp * ( ( ( a(:,ia+i) + a(:,ie+i) ) + a(:,ic+i) )                      &
+                                                - ( a(:,ib+i) + a(:,id+i) ) )
+             c(:,jc+j) = 2.0_wp * ( ( ( a(:,ia+i) + a(:,ie+i) ) - a(:,ic+i) )                      &
+                                                  - ( b(:,ib+i) - b(:,id+i) ) )
+             c(:,jg+j) = 2.0_wp * ( ( ( a(:,ia+i) + a(:,ie+i) ) - a(:,ic+i) )                      &
+                                                  + ( b(:,ib+i) - b(:,id+i) ) )
+             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) ) )
+             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) ) )
+             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) ) )
+             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) ) )
              ibase = ibase + inc1
              jbase = jbase + inc2
