Index: /palm/trunk/SOURCE/Makefile
===================================================================
--- /palm/trunk/SOURCE/Makefile	(revision 4365)
+++ /palm/trunk/SOURCE/Makefile	(revision 4366)
@@ -25,4 +25,7 @@
 # -----------------
 # $Id$
+# add dependency on fft for transpose
+# 
+# 4347 2019-12-18 13:18:33Z suehring
 # add dependency to basic_constants_and_equations_mod for dynamics_mod
 # 
@@ -1160,4 +1163,5 @@
 transpose.o: \
 	cpulog_mod.o \
+	fft_xy_mod.o \
 	mod_kinds.o \
 	modules.o
Index: /palm/trunk/SOURCE/fft_xy_mod.f90
===================================================================
--- /palm/trunk/SOURCE/fft_xy_mod.f90	(revision 4365)
+++ /palm/trunk/SOURCE/fft_xy_mod.f90	(revision 4366)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Vectorized Temperton-fft added
+! 
+! 4360 2020-01-07 11:25:50Z suehring
 ! Corrected "Former revisions" section
 ! 
@@ -41,4 +44,5 @@
 !> Fast Fourier transformation along x and y for 1d domain decomposition along x.
 !> Original version: Klaus Ketelsen (May 2002)
+!> @todo openmp support for vectorized Temperton fft
 !------------------------------------------------------------------------------!
  MODULE fft_xy
@@ -46,5 +50,5 @@
 
     USE control_parameters,                                                    &
-        ONLY:  fft_method, message_string
+        ONLY:  fft_method, loop_optimization, message_string
         
     USE cuda_fft_interfaces
@@ -52,5 +56,5 @@
     USE indices,                                                               &
         ONLY:  nx, ny, nz
-        
+
 #if defined( __cuda_fft )
     USE ISO_C_BINDING
@@ -72,10 +76,11 @@
 
     PRIVATE
-    PUBLIC fft_x, fft_x_1d, fft_y, fft_y_1d, fft_init, fft_x_m, fft_y_m
+    PUBLIC fft_x, fft_x_1d, fft_y, fft_y_1d, fft_init, fft_x_m, fft_y_m, f_vec, temperton_fft_vec
 
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE, SAVE ::  ifax_x  !<
     INTEGER(iwp), DIMENSION(:), ALLOCATABLE, SAVE ::  ifax_y  !<
 
-    LOGICAL, SAVE ::  init_fft = .FALSE.  !<
+    LOGICAL, SAVE ::  init_fft = .FALSE.           !<
+    LOGICAL, SAVE ::  temperton_fft_vec = .FALSE.  !<
 
     REAL(wp), SAVE ::  dnx      !<
@@ -86,4 +91,6 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  trigs_x  !< 
     REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  trigs_y  !<
+
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE, SAVE ::  f_vec
 
 #if defined( __ibm )
@@ -200,4 +207,11 @@
        ENDIF
 
+!
+!--    Switch to tell the Poisson-solver that the vectorized version of Temperton-fft is to be used.
+       IF ( fft_method == 'temperton-algorithm' .AND. loop_optimization == 'vector' )  THEN
+          temperton_fft_vec = .TRUE.
+       ENDIF
+
+
 #if defined( _OPENACC ) && defined( __cuda_fft )
        fft_method = 'system-specific'
@@ -270,4 +284,10 @@
           CALL set99( trigs_y, ifax_y, ny+1 )
 
+          IF ( temperton_fft_vec )  THEN
+             ALLOCATE( f_vec((nyn_x-nys_x+1)*(nzt_x-nzb_x+1),0:nx+2) )
+          ENDIF
+
+
+
        ELSEIF ( fft_method == 'fftw' )  THEN
 !
@@ -312,5 +332,5 @@
 !------------------------------------------------------------------------------!
  
-    SUBROUTINE fft_x( ar, direction, ar_2d )
+    SUBROUTINE fft_x( ar, direction, ar_2d, ar_inv )
 
 
@@ -325,4 +345,5 @@
        INTEGER(iwp) ::  j          !<
        INTEGER(iwp) ::  k          !<
+       INTEGER(iwp) ::  mm         !<
 
        LOGICAL ::  forward_fft !<
@@ -331,4 +352,10 @@
        REAL(wp), DIMENSION(nx+2)   ::  work1  !<
        
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE           ::  work_vec  !<
+       REAL(wp), DIMENSION(0:nx,nys_x:nyn_x), OPTIONAL ::  ar_2d     !<
+
+       REAL(wp), DIMENSION(nys_x:nyn_x,nzb_x:nzt_x,0:nx), OPTIONAL ::  ar_inv   !<
+       REAL(wp), DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x)           ::  ar       !<
+
 #if defined( __ibm )
        REAL(wp), DIMENSION(nau2) ::  aux2  !< 
@@ -337,13 +364,7 @@
        REAL(wp), DIMENSION(6*(nx+1)) ::  work2  !<
 #elif defined( __cuda_fft )
-       COMPLEX(dp), DIMENSION(0:(nx+1)/2,nys_x:nyn_x,nzb_x:nzt_x) ::           &
-          ar_tmp  !<
+       COMPLEX(dp), DIMENSION(0:(nx+1)/2,nys_x:nyn_x,nzb_x:nzt_x) ::  ar_tmp  !<
        !$ACC DECLARE CREATE(ar_tmp)
 #endif
-
-       REAL(wp), DIMENSION(0:nx,nys_x:nyn_x), OPTIONAL   ::                    &
-          ar_2d   !<
-       REAL(wp), DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x) ::                    &
-          ar      !<
 
 !
@@ -427,45 +448,133 @@
           IF ( forward_fft )  THEN
 
-             !$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
-             !$OMP DO
-             DO  k = nzb_x, nzt_x
-                DO  j = nys_x, nyn_x
-
-                   work(0:nx) = ar(0:nx,j,k)
-                   CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
-
-                   DO  i = 0, (nx+1)/2
-                      ar(i,j,k) = work(2*i)
-                   ENDDO
-                   DO  i = 1, (nx+1)/2 - 1
-                      ar(nx+1-i,j,k) = work(2*i+1)
-                   ENDDO
-
-                ENDDO
-             ENDDO
-             !$OMP END PARALLEL
-
-          ELSE
-
-             !$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
-             !$OMP DO
-             DO  k = nzb_x, nzt_x
-                DO  j = nys_x, nyn_x
-
-                   DO  i = 0, (nx+1)/2
-                      work(2*i) = ar(i,j,k)
-                   ENDDO
-                   DO  i = 1, (nx+1)/2 - 1
-                      work(2*i+1) = ar(nx+1-i,j,k)
-                   ENDDO
-                   work(1)    = 0.0_wp
-                   work(nx+2) = 0.0_wp
-
-                   CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
-                   ar(0:nx,j,k) = work(0:nx)
-
-                ENDDO
-             ENDDO
-             !$OMP END PARALLEL
+             IF ( .NOT. temperton_fft_vec )  THEN
+
+                !$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
+                !$OMP DO
+                DO  k = nzb_x, nzt_x
+                   DO  j = nys_x, nyn_x
+
+                      work(0:nx) = ar(0:nx,j,k)
+                      CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, -1 )
+
+                      DO  i = 0, (nx+1)/2
+                         ar(i,j,k) = work(2*i)
+                      ENDDO
+                      DO  i = 1, (nx+1)/2 - 1
+                         ar(nx+1-i,j,k) = work(2*i+1)
+                      ENDDO
+
+                   ENDDO
+                ENDDO
+                !$OMP END PARALLEL
+
+             ELSE
+
+!
+!--             Vector version of the Temperton-algorithm. Computes multiple 1-D FFT's.
+                ALLOCATE( work_vec( (nyn_x-nys_x+1)*(nzt_x-nzb_x+1),nx+2) )
+!
+!--             f_vec is already set in transpose_zx
+                CALL fft991cy_vec( f_vec, work_vec, trigs_x, ifax_x, nx+1, -1 )
+                DEALLOCATE( work_vec )
+
+                IF ( PRESENT( ar_inv ) )  THEN
+
+                   DO  k = nzb_x, nzt_x
+                      DO  j = nys_x, nyn_x
+                         mm = j-nys_x+1+(k-nzb_x)*(nyn_x-nys_x+1)
+                         DO  i = 0, (nx+1)/2
+                            ar_inv(j,k,i) = f_vec(mm,2*i)
+                         ENDDO
+                         DO  i = 1, (nx+1)/2-1
+                            ar_inv(j,k,nx+1-i) = f_vec(mm,2*i+1)
+                         ENDDO
+                      ENDDO
+                   ENDDO
+
+                ELSE
+
+                   DO  k = nzb_x, nzt_x
+                      DO  j = nys_x, nyn_x
+                         mm = j-nys_x+1+(k-nzb_x)*(nyn_x-nys_x+1)
+                         DO  i = 0, (nx+1)/2
+                            ar(i,j,k) = f_vec(mm,2*i)
+                         ENDDO
+                         DO  i = 1, (nx+1)/2-1
+                            ar(nx+1-i,j,k) = f_vec(mm,2*i+1)
+                         ENDDO
+                      ENDDO
+                   ENDDO
+
+                ENDIF
+
+             ENDIF
+
+          ELSE
+
+!
+!--          Backward fft
+             IF ( .NOT. temperton_fft_vec )  THEN
+
+                !$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
+                !$OMP DO
+                DO  k = nzb_x, nzt_x
+                   DO  j = nys_x, nyn_x
+
+                      DO  i = 0, (nx+1)/2
+                         work(2*i) = ar(i,j,k)
+                      ENDDO
+                      DO  i = 1, (nx+1)/2 - 1
+                         work(2*i+1) = ar(nx+1-i,j,k)
+                      ENDDO
+                      work(1)    = 0.0_wp
+                      work(nx+2) = 0.0_wp
+
+                      CALL fft991cy( work, work1, trigs_x, ifax_x, 1, nx+1, nx+1, 1, 1 )
+                      ar(0:nx,j,k) = work(0:nx)
+
+                   ENDDO
+                ENDDO
+                !$OMP END PARALLEL
+
+             ELSE
+
+                IF ( PRESENT( ar_inv ) )  THEN
+
+                   DO  k = nzb_x, nzt_x
+                      DO  j = nys_x, nyn_x
+                         mm = j-nys_x+1+(k-nzb_x)*(nyn_x-nys_x+1)
+                         DO  i = 0, (nx+1)/2
+                            f_vec(mm,2*i) = ar_inv(j,k,i)
+                         ENDDO
+                         DO  i = 1, (nx+1)/2-1
+                            f_vec(mm,2*i+1) = ar_inv(j,k,nx+1-i)
+                         ENDDO
+                      ENDDO
+                   ENDDO
+
+                ELSE
+
+                   DO  k = nzb_x, nzt_x
+                      DO  j = nys_x, nyn_x
+                         mm = j-nys_x+1+(k-nzb_x)*(nyn_x-nys_x+1)
+                         DO  i = 0, (nx+1)/2
+                            f_vec(mm,2*i) = ar(i,j,k)
+                         ENDDO
+                         DO  i = 1, (nx+1)/2-1
+                            f_vec(mm,2*i+1) = ar(nx+1-i,j,k)
+                         ENDDO
+                      ENDDO
+                   ENDDO
+
+                ENDIF
+                f_vec(:,1)    = 0.0_wp
+                f_vec(:,nx+2) = 0.0_wp
+
+                ALLOCATE( work_vec((nyn_x-nys_x+1)*(nzt_x-nzb_x+1),nx+2) )
+                CALL fft991cy_vec( f_vec, work_vec, trigs_x, ifax_x, nx+1, 1 )
+                DEALLOCATE( work_vec )
+
+             ENDIF
 
           ENDIF
@@ -941,5 +1050,5 @@
  
     SUBROUTINE fft_y( ar, direction, ar_tr, nxl_y_bound, nxr_y_bound, nxl_y_l, &
-                      nxr_y_l )
+                      nxr_y_l, ar_inv )
 
 
@@ -952,4 +1061,5 @@
        INTEGER(iwp) ::  jshape(1)    !<
        INTEGER(iwp) ::  k            !<
+       INTEGER(iwp) ::  mm           !<
        INTEGER(iwp) ::  nxl_y_bound  !<
        INTEGER(iwp) ::  nxl_y_l      !<
@@ -962,4 +1072,11 @@
        REAL(wp), DIMENSION(ny+2)   ::  work1  !<
        
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  f_vec
+       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  work_vec
+
+       REAL(wp), DIMENSION(0:ny,nxl_y_l:nxr_y_l,nzb_y:nzt_y)                   ::  ar      !<
+       REAL(wp), DIMENSION(nxl_y:nxr_y,nzb_y:nzt_y,0:ny), OPTIONAL             ::  ar_inv  !<
+       REAL(wp), DIMENSION(0:ny,nxl_y_bound:nxr_y_bound,nzb_y:nzt_y), OPTIONAL ::  ar_tr   !<
+
        COMPLEX(wp), DIMENSION(:), ALLOCATABLE ::  cwork  !<
        
@@ -975,8 +1092,4 @@
 #endif
 
-       REAL(wp), DIMENSION(0:ny,nxl_y_l:nxr_y_l,nzb_y:nzt_y)         ::        &
-          ar     !<
-       REAL(wp), DIMENSION(0:ny,nxl_y_bound:nxr_y_bound,nzb_y:nzt_y) ::        &
-          ar_tr  !<
 
        IF ( direction == 'forward' )  THEN
@@ -1057,45 +1170,153 @@
           IF ( forward_fft )  THEN
 
-             !$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
-             !$OMP DO
-             DO  k = nzb_y, nzt_y
-                DO  i = nxl_y_l, nxr_y_l
-
-                   work(0:ny) = ar(0:ny,i,k)
-                   CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
-
-                   DO  j = 0, (ny+1)/2
-                      ar_tr(j,i,k) = work(2*j)
-                   ENDDO
-                   DO  j = 1, (ny+1)/2 - 1
-                      ar_tr(ny+1-j,i,k) = work(2*j+1)
-                   ENDDO
-
-                ENDDO
-             ENDDO
-             !$OMP END PARALLEL
-
-          ELSE
-
-             !$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
-             !$OMP DO
-             DO  k = nzb_y, nzt_y
-                DO  i = nxl_y_l, nxr_y_l
-
-                   DO  j = 0, (ny+1)/2
-                      work(2*j) = ar_tr(j,i,k)
-                   ENDDO
-                   DO  j = 1, (ny+1)/2 - 1
-                      work(2*j+1) = ar_tr(ny+1-j,i,k)
-                   ENDDO
-                   work(1)    = 0.0_wp
-                   work(ny+2) = 0.0_wp
-
-                   CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
-                   ar(0:ny,i,k) = work(0:ny)
-
-                ENDDO
-             ENDDO
-             !$OMP END PARALLEL
+             IF ( .NOT. temperton_fft_vec )  THEN
+
+                !$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
+                !$OMP DO
+                DO  k = nzb_y, nzt_y
+                   DO  i = nxl_y_l, nxr_y_l
+
+                      work(0:ny) = ar(0:ny,i,k)
+                      CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, -1 )
+
+                      DO  j = 0, (ny+1)/2
+                         ar_tr(j,i,k) = work(2*j)
+                      ENDDO
+                      DO  j = 1, (ny+1)/2 - 1
+                         ar_tr(ny+1-j,i,k) = work(2*j+1)
+                      ENDDO
+
+                   ENDDO
+                ENDDO
+                !$OMP END PARALLEL
+
+             ELSE
+!
+!--             Vector version of Temperton-fft. Computes multiple 1-D FFT's.
+                ALLOCATE( f_vec((nxr_y_l-nxl_y_l+1)*(nzt_y-nzb_y+1),0:ny+2) )
+
+                mm = 1
+                DO  k = nzb_y, nzt_y
+                   DO  i = nxl_y_l, nxr_y_l
+                      f_vec(mm,0:nx) = ar(0:nx,i,k)
+                      mm = mm+1
+                   ENDDO
+                ENDDO
+
+                ALLOCATE( work_vec( (nxr_y_l-nxl_y_l+1)*(nzt_y-nzb_y+1),ny+2) )
+                CALL fft991cy_vec( f_vec, work_vec, trigs_y, ifax_y, ny+1, -1 )
+                DEALLOCATE( work_vec )
+
+                IF( PRESENT( ar_inv ) )  THEN
+
+                   DO  k = nzb_y, nzt_y
+                      DO  i = nxl_y_l, nxr_y_l
+                         mm = i-nxl_y_l+1+(k-nzb_y)*(nxr_y_l-nxl_y_l+1)
+                         DO  j = 0, (ny+1)/2
+                            ar_inv(i,k,j) = f_vec(mm,2*j)
+                         ENDDO
+                         DO  j = 1, (ny+1)/2 - 1
+                            ar_inv(i,k,ny+1-j) = f_vec(mm,2*j+1)
+                         ENDDO
+                      ENDDO
+                   ENDDO
+
+                ELSE
+
+                   DO  k = nzb_y, nzt_y
+                      DO  i = nxl_y_l, nxr_y_l
+                         mm = i-nxl_y_l+1+(k-nzb_y)*(nxr_y_l-nxl_y_l+1)
+                         DO  j = 0, (ny+1)/2
+                            ar(j,i,k) = f_vec(mm,2*j)
+                         ENDDO
+                         DO  j = 1, (ny+1)/2 - 1
+                            ar(ny+1-j,i,k) = f_vec(mm,2*j+1)
+                         ENDDO
+                      ENDDO
+                   ENDDO
+
+                ENDIF
+
+                DEALLOCATE( f_vec )
+
+             ENDIF
+
+          ELSE
+
+             IF ( .NOT. temperton_fft_vec )  THEN
+
+                !$OMP PARALLEL PRIVATE ( work, work1, i, j, k )
+                !$OMP DO
+                DO  k = nzb_y, nzt_y
+                   DO  i = nxl_y_l, nxr_y_l
+
+                      DO  j = 0, (ny+1)/2
+                         work(2*j) = ar_tr(j,i,k)
+                      ENDDO
+                      DO  j = 1, (ny+1)/2 - 1
+                         work(2*j+1) = ar_tr(ny+1-j,i,k)
+                      ENDDO
+                      work(1)    = 0.0_wp
+                      work(ny+2) = 0.0_wp
+
+                      CALL fft991cy( work, work1, trigs_y, ifax_y, 1, ny+1, ny+1, 1, 1 )
+                      ar(0:ny,i,k) = work(0:ny)
+
+                   ENDDO
+                ENDDO
+                !$OMP END PARALLEL
+
+             ELSE
+
+                ALLOCATE( f_vec((nxr_y_l-nxl_y_l+1)*(nzt_y-nzb_y+1),0:ny+2) )
+
+                IF ( PRESENT( ar_inv ) )  THEN
+
+                   DO  k = nzb_y, nzt_y
+                      DO  i = nxl_y_l, nxr_y_l
+                         mm = i-nxl_y_l+1+(k-nzb_y)*(nxr_y_l-nxl_y_l+1)
+                         DO  j = 0, (ny+1)/2
+                            f_vec(mm,2*j) = ar_inv(i,k,j)
+                         ENDDO
+                         DO  j = 1, (ny+1)/2 - 1
+                            f_vec(mm,2*j+1) = ar_inv(i,k,ny+1-j)
+                         ENDDO
+                      ENDDO
+                   ENDDO
+
+                ELSE
+
+                   DO  k = nzb_y, nzt_y
+                      DO  i = nxl_y_l, nxr_y_l
+                         mm = i-nxl_y_l+1+(k-nzb_y)*(nxr_y_l-nxl_y_l+1)
+                         DO  j = 0, (ny+1)/2
+                            f_vec(mm,2*j) = ar(j,i,k)
+                         ENDDO
+                         DO  j = 1, (ny+1)/2 - 1
+                            f_vec(mm,2*j+1) = ar(ny+1-j,i,k)
+                         ENDDO
+                      ENDDO
+                   ENDDO
+
+                ENDIF
+
+                f_vec(:,1)    = 0.0_wp
+                f_vec(:,ny+2) = 0.0_wp
+
+                ALLOCATE( work_vec((nxr_y_l-nxl_y_l+1)*(nzt_y-nzb_y+1),ny+2) )
+                CALL fft991cy_vec( f_vec, work_vec, trigs_y, ifax_y, ny+1, 1 )
+                DEALLOCATE( work_vec )
+
+                mm = 1
+                DO  k = nzb_y, nzt_y
+                   DO  i = nxl_y_l, nxr_y_l
+                      ar(0:ny,i,k) = f_vec(mm,0:ny)
+                      mm = mm+1
+                   ENDDO
+                ENDDO
+
+                DEALLOCATE( f_vec )
+
+             ENDIF
 
           ENDIF
Index: /palm/trunk/SOURCE/poisfft_mod.f90
===================================================================
--- /palm/trunk/SOURCE/poisfft_mod.f90	(revision 4365)
+++ /palm/trunk/SOURCE/poisfft_mod.f90	(revision 4366)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! modification concerning NEC vectorizatio
+! 
+! 4360 2020-01-07 11:25:50Z suehring
 ! Corrected "Former revisions" section
 ! 
@@ -50,5 +53,6 @@
 
     USE fft_xy,                                                                &
-        ONLY:  fft_init, fft_y, fft_y_1d, fft_y_m, fft_x, fft_x_1d, fft_x_m
+        ONLY:  fft_init, fft_y, fft_y_1d, fft_y_m, fft_x, fft_x_1d, fft_x_m,   &
+               temperton_fft_vec
 
     USE indices,                                                               &
@@ -207,5 +211,12 @@
 
           CALL cpu_log( log_point_s(4), 'fft_x', 'start' )
-          CALL fft_x( ar, 'forward' )
+          IF ( temperton_fft_vec )  THEN
+!
+!--          Vector version outputs a transformed array ar_inv that does not require resorting
+!--          (which is done for ar further below)
+             CALL fft_x( ar, 'forward',  ar_inv=ar_inv)
+          ELSE
+             CALL fft_x( ar, 'forward')
+          ENDIF
           CALL cpu_log( log_point_s(4), 'fft_x', 'pause' )
 
@@ -213,12 +224,19 @@
 !--       Transposition x --> y
           CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
-          CALL resort_for_xy( ar, ar_inv )
+          IF( .NOT. temperton_fft_vec )  CALL resort_for_xy( ar, ar_inv )
           CALL transpose_xy( ar_inv, ar )
           CALL cpu_log( log_point_s(5), 'transpo forward', 'pause' )
 
           CALL cpu_log( log_point_s(7), 'fft_y', 'start' )
-          CALL fft_y( ar, 'forward', ar_tr = ar,                &
-                      nxl_y_bound = nxl_y, nxr_y_bound = nxr_y, &
-                      nxl_y_l = nxl_y, nxr_y_l = nxr_y )
+          IF ( temperton_fft_vec )  THEN
+!
+!--          Input array ar_inv from fft_x can be directly used here.
+!--          The output (also in array ar_inv) does not require resorting below.
+             CALL fft_y( ar, 'forward', ar_inv = ar_inv, nxl_y_bound = nxl_y, nxr_y_bound = nxr_y, &
+                         nxl_y_l = nxl_y, nxr_y_l = nxr_y )
+          ELSE
+             CALL fft_y( ar, 'forward', ar_tr = ar, nxl_y_bound = nxl_y, nxr_y_bound = nxr_y,      &
+                         nxl_y_l = nxl_y, nxr_y_l = nxr_y )
+          ENDIF
           CALL cpu_log( log_point_s(7), 'fft_y', 'pause' )
 
@@ -226,5 +244,5 @@
 !--       Transposition y --> z
           CALL cpu_log( log_point_s(5), 'transpo forward', 'continue' )
-          CALL resort_for_yz( ar, ar_inv )
+          IF ( .NOT. temperton_fft_vec )  CALL resort_for_yz( ar, ar_inv )
           CALL transpose_yz( ar_inv, ar )
           CALL cpu_log( log_point_s(5), 'transpo forward', 'stop' )
@@ -241,11 +259,21 @@
           CALL cpu_log( log_point_s(8), 'transpo invers', 'start' )
           CALL transpose_zy( ar, ar_inv )
-          CALL resort_for_zy( ar_inv, ar )
+!
+!--       The fft_y below (vector branch) can directly process ar_inv (i.e. does not require a
+!--       resorting)
+          IF ( .NOT. temperton_fft_vec )  CALL resort_for_zy( ar_inv, ar )
           CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
 
           CALL cpu_log( log_point_s(7), 'fft_y', 'continue' )
-          CALL fft_y( ar, 'backward', ar_tr = ar,               &
-                      nxl_y_bound = nxl_y, nxr_y_bound = nxr_y, &
-                      nxl_y_l = nxl_y, nxr_y_l = nxr_y )
+          IF ( temperton_fft_vec )  THEN
+!
+!--          Output array ar_inv can be used as input to the below fft_x routine without resorting
+             CALL fft_y( ar, 'backward', ar_inv = ar_inv, nxl_y_bound = nxl_y, nxr_y_bound = nxr_y,&
+                         nxl_y_l = nxl_y, nxr_y_l = nxr_y )
+          ELSE
+             CALL fft_y( ar, 'backward', ar_tr = ar, nxl_y_bound = nxl_y, nxr_y_bound = nxr_y,     &
+                         nxl_y_l = nxl_y, nxr_y_l = nxr_y )
+          ENDIF
+
           CALL cpu_log( log_point_s(7), 'fft_y', 'stop' )
 
@@ -254,9 +282,13 @@
           CALL cpu_log( log_point_s(8), 'transpo invers', 'continue' )
           CALL transpose_yx( ar, ar_inv )
-          CALL resort_for_yx( ar_inv, ar )
+          IF ( .NOT. temperton_fft_vec )  CALL resort_for_yx( ar_inv, ar )
           CALL cpu_log( log_point_s(8), 'transpo invers', 'pause' )
 
           CALL cpu_log( log_point_s(4), 'fft_x', 'continue' )
-          CALL fft_x( ar, 'backward' )
+          IF ( temperton_fft_vec )  THEN
+             CALL fft_x( ar, 'backward',  ar_inv=ar_inv )
+          ELSE
+             CALL fft_x( ar, 'backward' )
+          ENDIF
           CALL cpu_log( log_point_s(4), 'fft_x', 'stop' )
 
Index: /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90
===================================================================
--- /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 4365)
+++ /palm/trunk/SOURCE/surface_layer_fluxes_mod.f90	(revision 4366)
@@ -26,4 +26,7 @@
 ! -----------------
 ! $Id$
+! vector version for calculation of Obukhov length via Newton iteration added
+! 
+! 4360 2020-01-07 11:25:50Z suehring
 ! Calculation of diagnostic-only 2-m potential temperature moved to 
 ! diagnostic_output_quantities.
@@ -112,5 +115,5 @@
                constant_waterflux, coupling_mode,                              &
                debug_output_timestep,                                          &
-               humidity,                                                       &
+               humidity, loop_optimization,                                    &
                ibc_e_b, ibc_pt_b, indoor_model,                                &
                land_surface, large_scale_forcing, lsf_surf, message_string,    &
@@ -859,4 +862,6 @@
        INTEGER(iwp) ::  m       !< loop variable over all horizontal wall elements 
 
+       LOGICAL, DIMENSION(surf%ns) ::  convergence_reached  !< convergence switch for vectorization
+
        REAL(wp)     :: f,      & !< Function for Newton iteration: f = Ri - [...]/[...]^2 = 0
                        f_d_ol, & !< Derivative of f
@@ -866,4 +871,7 @@
                        ol_u      !< Upper bound of L for Newton iteration
 
+       REAL(wp), DIMENSION(surf%ns) ::  ol_old_vec  !< temporary array required for vectorization
+       REAL(wp), DIMENSION(surf%ns) ::  z_mo_vec    !< temporary array required for vectorization
+
 !
 !--    Evaluate bulk Richardson number (calculation depends on
@@ -926,134 +934,265 @@
        ENDIF
 
-!
-!--    Calculate the Obukhov length using Newton iteration
-       !$OMP PARALLEL DO PRIVATE(i, j, z_mo) &
-       !$OMP PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol)
-       !$ACC PARALLEL LOOP PRIVATE(i, j, z_mo) &
-       !$ACC PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) &
-       !$ACC PRESENT(surf)
-       DO  m = 1, surf%ns
-
-          i   = surf%i(m)            
-          j   = surf%j(m)
-
-          z_mo = surf%z_mo(m)
-
-!
-!--       Store current value in case the Newton iteration fails
-          ol_old = surf%ol(m)
-
-!
-!--       Ensure that the bulk Richardson number and the Obukhov 
-!--       length have the same sign
-          IF ( surf%rib(m) * surf%ol(m) < 0.0_wp  .OR.                      &
-               ABS( surf%ol(m) ) == ol_max )  THEN
-             IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) =  0.01_wp
-             IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp
-          ENDIF
+       IF ( loop_optimization == 'cache' )  THEN
+!
+!--       Calculate the Obukhov length using Newton iteration
+          !$OMP PARALLEL DO PRIVATE(i, j, z_mo) &
+          !$OMP PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol)
+          !$ACC PARALLEL LOOP PRIVATE(i, j, z_mo) &
+          !$ACC PRIVATE(ol_old, ol_m, ol_l, ol_u, f, f_d_ol) &
+          !$ACC PRESENT(surf)
+          DO  m = 1, surf%ns
+
+             i   = surf%i(m)
+             j   = surf%j(m)
+
+             z_mo = surf%z_mo(m)
+
+!
+!--          Store current value in case the Newton iteration fails
+             ol_old = surf%ol(m)
+
+!
+!--          Ensure that the bulk Richardson number and the Obukhov
+!--          length have the same sign
+             IF ( surf%rib(m) * surf%ol(m) < 0.0_wp  .OR.  ABS( surf%ol(m) ) == ol_max )  THEN
+                IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) =  0.01_wp
+                IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp
+             ENDIF
+!
+!--          Iteration to find Obukhov length
+             iter = 0
+             DO
+                iter = iter + 1
+!
+!--             In case of divergence, use the value of the previous time step
+                IF ( iter > 1000 )  THEN
+                   surf%ol(m) = ol_old
+                   EXIT
+                ENDIF
+
+                ol_m = surf%ol(m)
+                ol_l = ol_m - 0.001_wp * ol_m
+                ol_u = ol_m + 0.001_wp * ol_m
+
+
+                IF ( ibc_pt_b /= 1 )  THEN
+!
+!--                Calculate f = Ri - [...]/[...]^2 = 0
+                   f = surf%rib(m) - ( z_mo / ol_m ) * ( LOG( z_mo / surf%z0h(m) )                 &
+                                                       - psi_h( z_mo / ol_m )                      &
+                                                       + psi_h( surf%z0h(m) / ol_m )               &
+                                                       ) /                                         &
+                                                       ( LOG( z_mo / surf%z0(m) )                  &
+                                                      - psi_m( z_mo / ol_m )                       &
+                                                      + psi_m( surf%z0(m) /  ol_m )                &
+                                                       )**2
+
+!
+!--                Calculate df/dL
+                   f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo / surf%z0h(m) )                      &
+                                                  - psi_h( z_mo / ol_u )                           &
+                                                  + psi_h( surf%z0h(m) / ol_u )                    &
+                                                  ) /                                              &
+                                                  ( LOG( z_mo / surf%z0(m) )                       &
+                                                  - psi_m( z_mo / ol_u )                           &
+                                                  + psi_m( surf%z0(m) / ol_u )                     &
+                                                  )**2                                             &
+                              + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) )                      &
+                                                  - psi_h( z_mo / ol_l )                           &
+                                                  + psi_h( surf%z0h(m) / ol_l )                    &
+                                                  ) /                                              &
+                                                  ( LOG( z_mo / surf%z0(m) )                       &
+                                                  - psi_m( z_mo / ol_l )                           &
+                                                  + psi_m( surf%z0(m) / ol_l )                     &
+                                                  )**2                                             &
+                           ) / ( ol_u - ol_l )
+                ELSE
+!
+!--                Calculate f = Ri - 1 /[...]^3 = 0
+                   f = surf%rib(m) - ( z_mo / ol_m ) / ( LOG( z_mo / surf%z0(m) )                  &
+                                                       - psi_m( z_mo / ol_m )                      &
+                                                       + psi_m( surf%z0(m) / ol_m )                &
+                                                       )**3
+
+!
+!--                Calculate df/dL
+                   f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo / surf%z0(m) )                       &
+                                                  - psi_m( z_mo / ol_u )                           &
+                                                  + psi_m( surf%z0(m) / ol_u )                     &
+                                                  )**3                                             &
+                              + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) )                       &
+                                                  - psi_m( z_mo / ol_l )                           &
+                                                  + psi_m( surf%z0(m) / ol_l )                     &
+                                                  )**3                                             &
+                             ) / ( ol_u - ol_l )
+                ENDIF
+!
+!--             Calculate new L
+                surf%ol(m) = ol_m - f / f_d_ol
+
+!
+!--             Ensure that the bulk Richardson number and the Obukhov
+!--             length have the same sign and ensure convergence.
+                IF ( surf%ol(m) * ol_m < 0.0_wp )  surf%ol(m) = ol_m * 0.5_wp
+!
+!--             If unrealistic value occurs, set L to the maximum
+!--             value that is allowed
+                IF ( ABS( surf%ol(m) ) > ol_max )  THEN
+                   surf%ol(m) = ol_max
+                   EXIT
+                ENDIF
+!
+!--             Assure that Obukhov length does not become zero. If the limit is
+!--             reached, exit the loop.
+                IF ( ABS( surf%ol(m) ) < 1E-5_wp )  THEN
+                   surf%ol(m) = SIGN( 1E-5_wp, surf%ol(m) )
+                   EXIT
+                ENDIF
+!
+!--             Check for convergence
+                IF ( ABS( ( surf%ol(m) - ol_m ) /  surf%ol(m) ) < 1.0E-4_wp )  EXIT
+
+             ENDDO
+          ENDDO
+
+!
+!--    Vector Version
+       ELSE
+!
+!--       Calculate the Obukhov length using Newton iteration
+!--       First set arrays required for vectorization
+          DO  m = 1, surf%ns
+
+             z_mo_vec(m) = surf%z_mo(m)
+
+!
+!--          Store current value in case the Newton iteration fails
+             ol_old_vec(m) = surf%ol(m)
+
+!
+!--          Ensure that the bulk Richardson number and the Obukhov length have the same sign
+             IF ( surf%rib(m) * surf%ol(m) < 0.0_wp  .OR.  ABS( surf%ol(m) ) == ol_max )  THEN
+                IF ( surf%rib(m) > 1.0_wp ) surf%ol(m) =  0.01_wp
+                IF ( surf%rib(m) < 0.0_wp ) surf%ol(m) = -0.01_wp
+             ENDIF
+
+          ENDDO
+
 !
 !--       Iteration to find Obukhov length
+          convergence_reached(:) = .FALSE.
           iter = 0
           DO
+
              iter = iter + 1
 !
-!--          In case of divergence, use the value of the previous time step
+!--          In case of divergence, use the value(s) of the previous time step
              IF ( iter > 1000 )  THEN
-                surf%ol(m) = ol_old
+                DO  m = 1, surf%ns
+                   IF ( .NOT. convergence_reached(m) )  surf%ol(1:surf%ns) = ol_old
+                ENDDO
                 EXIT
              ENDIF
 
-             ol_m = surf%ol(m)
-             ol_l = ol_m - 0.001_wp * ol_m
-             ol_u = ol_m + 0.001_wp * ol_m
-
-
-             IF ( ibc_pt_b /= 1 )  THEN
-!
-!--             Calculate f = Ri - [...]/[...]^2 = 0
-                f = surf%rib(m) - ( z_mo / ol_m ) * (                          &
-                                              LOG( z_mo / surf%z0h(m) )        &
-                                              - psi_h( z_mo / ol_m )           &
-                                              + psi_h( surf%z0h(m) /           &
-                                                       ol_m )                  &  
-                                                     )                         &
-                                           / ( LOG( z_mo / surf%z0(m) )        &
-                                              - psi_m( z_mo / ol_m )           &
-                                              + psi_m( surf%z0(m) /  ol_m )    &
-                                                 )**2
-
-!
-!--              Calculate df/dL
-                 f_d_ol = ( - ( z_mo / ol_u ) * ( LOG( z_mo /                  &
-                                                          surf%z0h(m) )        &
-                                         - psi_h( z_mo / ol_u )                &
-                                         + psi_h( surf%z0h(m) / ol_u )         &
-                                           )                                   &
-                                         / ( LOG( z_mo / surf%z0(m) )          &
-                                         - psi_m( z_mo / ol_u )                &
-                                         + psi_m( surf%z0(m) / ol_u )          &
-                                           )**2                                &
-                        + ( z_mo / ol_l ) * ( LOG( z_mo / surf%z0h(m) )        &
-                                         - psi_h( z_mo / ol_l )                &
-                                         + psi_h( surf%z0h(m) / ol_l )         &
-                                           )                                   &
-                                         / ( LOG( z_mo / surf%z0(m) )          &
-                                         - psi_m( z_mo / ol_l )                &
-                                         + psi_m( surf%z0(m) / ol_l )          &
-                                           )**2                                &
-                          ) / ( ol_u - ol_l )
-             ELSE
-!
-!--             Calculate f = Ri - 1 /[...]^3 = 0
-                f = surf%rib(m) - ( z_mo / ol_m ) /                            &
-                                             ( LOG( z_mo / surf%z0(m) )        &
-                                         - psi_m( z_mo / ol_m )                &
-                                         + psi_m( surf%z0(m) / ol_m )          &
-                                             )**3
-
-!
-!--             Calculate df/dL
-                f_d_ol = ( - ( z_mo / ol_u ) / ( LOG( z_mo / surf%z0(m) )      &
-                                         - psi_m( z_mo / ol_u )                &
-                                         + psi_m( surf%z0(m) / ol_u )          &
-                                                  )**3                         &
-                           + ( z_mo / ol_l ) / ( LOG( z_mo / surf%z0(m) )      &
-                                         - psi_m( z_mo / ol_l )                &
-                                         + psi_m( surf%z0(m) / ol_l )          &
-                                            )**3                               &
-                          ) / ( ol_u - ol_l )
-             ENDIF
-!
-!--          Calculate new L
-             surf%ol(m) = ol_m - f / f_d_ol
-
-!
-!--          Ensure that the bulk Richardson number and the Obukhov 
-!--          length have the same sign and ensure convergence.
-             IF ( surf%ol(m) * ol_m < 0.0_wp )  surf%ol(m) = ol_m * 0.5_wp
-!
-!--          If unrealistic value occurs, set L to the maximum
-!--          value that is allowed
-             IF ( ABS( surf%ol(m) ) > ol_max )  THEN
-                surf%ol(m) = ol_max
-                EXIT
-             ENDIF
-!
-!--          Assure that Obukhov length does not become zero. If the limit is
-!--          reached, exit the loop.
-             IF ( ABS( surf%ol(m) ) < 1E-5_wp )  THEN
-                surf%ol(m) = SIGN( 1E-5_wp, surf%ol(m) )
-                EXIT
-             ENDIF
-!
-!--          Check for convergence
-             IF ( ABS( ( surf%ol(m) - ol_m ) /  surf%ol(m) ) < 1.0E-4_wp )  THEN
-                EXIT
-             ELSE
-                CYCLE
-             ENDIF              
-
-          ENDDO
-       ENDDO
+
+             DO  m = 1, surf%ns
+
+                IF ( convergence_reached(m) )  CYCLE
+
+                ol_m = surf%ol(m)
+                ol_l = ol_m - 0.001_wp * ol_m
+                ol_u = ol_m + 0.001_wp * ol_m
+
+
+                IF ( ibc_pt_b /= 1 )  THEN
+!
+!--                Calculate f = Ri - [...]/[...]^2 = 0
+                   f = surf%rib(m) - ( z_mo_vec(m) / ol_m ) * ( LOG( z_mo_vec(m) / surf%z0h(m) )   &
+                                                              - psi_h( z_mo_vec(m) / ol_m )        &
+                                                              + psi_h( surf%z0h(m) / ol_m )        &
+                                                              ) /                                  &
+                                                              ( LOG( z_mo_vec(m) / surf%z0(m) )    &
+                                                             - psi_m( z_mo_vec(m) / ol_m )         &
+                                                             + psi_m( surf%z0(m) /  ol_m )         &
+                                                              )**2
+
+!
+!--                Calculate df/dL
+                   f_d_ol = ( - ( z_mo_vec(m) / ol_u ) * ( LOG( z_mo_vec(m) / surf%z0h(m) )        &
+                                                         - psi_h( z_mo_vec(m) / ol_u )             &
+                                                         + psi_h( surf%z0h(m) / ol_u )             &
+                                                         ) /                                       &
+                                                         ( LOG( z_mo_vec(m) / surf%z0(m) )         &
+                                                         - psi_m( z_mo_vec(m) / ol_u )             &
+                                                         + psi_m( surf%z0(m) / ol_u )              &
+                                                         )**2                                      &
+                              + ( z_mo_vec(m) / ol_l ) * ( LOG( z_mo_vec(m) / surf%z0h(m) )        &
+                                                         - psi_h( z_mo_vec(m) / ol_l )             &
+                                                         + psi_h( surf%z0h(m) / ol_l )             &
+                                                         ) /                                       &
+                                                         ( LOG( z_mo_vec(m) / surf%z0(m) )         &
+                                                         - psi_m( z_mo_vec(m) / ol_l )             &
+                                                         + psi_m( surf%z0(m) / ol_l )              &
+                                                         )**2                                      &
+                            ) / ( ol_u - ol_l )
+                ELSE
+!
+!--                Calculate f = Ri - 1 /[...]^3 = 0
+                   f = surf%rib(m) - ( z_mo_vec(m) / ol_m ) / ( LOG( z_mo_vec(m) / surf%z0(m) )    &
+                                                              - psi_m( z_mo_vec(m) / ol_m )        &
+                                                              + psi_m( surf%z0(m) / ol_m )         &
+                                                              )**3
+
+!
+!--                Calculate df/dL
+                   f_d_ol = ( - ( z_mo_vec(m) / ol_u ) / ( LOG( z_mo_vec(m) / surf%z0(m) )         &
+                                                         - psi_m( z_mo_vec(m) / ol_u )             &
+                                                         + psi_m( surf%z0(m) / ol_u )              &
+                                                         )**3                                      &
+                              + ( z_mo_vec(m) / ol_l ) / ( LOG( z_mo_vec(m) / surf%z0(m) )         &
+                                                         - psi_m( z_mo_vec(m) / ol_l )             &
+                                                         + psi_m( surf%z0(m) / ol_l )              &
+                                                         )**3                                      &
+                            ) / ( ol_u - ol_l )
+                ENDIF
+!
+!--             Calculate new L
+                surf%ol(m) = ol_m - f / f_d_ol
+
+!
+!--             Ensure that the bulk Richardson number and the Obukhov
+!--             length have the same sign and ensure convergence.
+                IF ( surf%ol(m) * ol_m < 0.0_wp )  surf%ol(m) = ol_m * 0.5_wp
+
+!
+!--             Check for convergence
+!--             This check does not modify surf%ol, therefore this is done first
+                IF ( ABS( ( surf%ol(m) - ol_m ) /  surf%ol(m) ) < 1.0E-4_wp )  THEN
+                   convergence_reached(m) = .TRUE.
+                ENDIF
+!
+!--             If unrealistic value occurs, set L to the maximum allowed value
+                IF ( ABS( surf%ol(m) ) > ol_max )  THEN
+                   surf%ol(m) = ol_max
+                   convergence_reached(m) = .TRUE.
+                ENDIF
+
+             ENDDO
+!
+!--          Assure that Obukhov length does not become zero
+             DO  m = 1, surf%ns
+                IF ( convergence_reached(m) )  CYCLE
+                IF ( ABS( surf%ol(m) ) < 1E-5_wp )  THEN
+                   surf%ol(m) = SIGN( 10E-6_wp, surf%ol(m) )
+                   convergence_reached(m) = .TRUE.
+                ENDIF
+             ENDDO
+
+             IF ( ALL( convergence_reached ) )  EXIT
+
+          ENDDO  ! end of iteration loop
+
+       ENDIF  ! end of vector branch
 
     END SUBROUTINE calc_ol
Index: /palm/trunk/SOURCE/surface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/surface_mod.f90	(revision 4365)
+++ /palm/trunk/SOURCE/surface_mod.f90	(revision 4366)
@@ -26,4 +26,7 @@
 ! -----------------
 ! $Id$
+! workaround implemented to avoid vectorization bug on NEC Aurora
+! 
+! 4360 2020-01-07 11:25:50Z suehring
 ! Fix also remaining message calls.
 ! 
@@ -2375,4 +2378,11 @@
                                                  num_def_h_kji(0) - 1
              start_index_def_h(0)           = surf_def_h(0)%end_index(j,i) + 1
+!
+!--          ATTENTION:
+!--          workaround to prevent vectorization bug on NEC Aurora
+             IF ( start_index_def_h(0) < -99999 )  THEN
+                PRINT*, 'i=', i, ' j=',j, ' s=',surf_def_h(0)%start_index(j,i),                    &
+                        ' e=', surf_def_h(0)%end_index(j,i)
+             ENDIF
 !
 !--          Downward-facing surfaces, except model top
Index: /palm/trunk/SOURCE/temperton_fft_mod.f90
===================================================================
--- /palm/trunk/SOURCE/temperton_fft_mod.f90	(revision 4365)
+++ /palm/trunk/SOURCE/temperton_fft_mod.f90	(revision 4366)
@@ -9,4 +9,7 @@
 ! -----------------
 ! $Id$
+! vectorized routines added
+! 
+! 4182 2019-08-22 15:20:23Z scharf
 ! Corrected "Former revisions" section
 ! 
@@ -30,5 +33,19 @@
     PRIVATE
 
-    PUBLIC set99, fft991cy
+!
+!-- No interfaces for the serial routines, because these are still writte in FORTRAN77
+    INTERFACE fft991cy_vec
+       MODULE PROCEDURE fft991cy_vec
+    END INTERFACE fft991cy_vec
+
+    INTERFACE qpassm_vec
+       MODULE PROCEDURE qpassm_vec
+    END INTERFACE qpassm_vec
+
+    INTERFACE rpassm_vec
+       MODULE PROCEDURE rpassm_vec
+    END INTERFACE rpassm_vec
+
+    PUBLIC set99, fft991cy, fft991cy_vec
 
     INTEGER(iwp), PARAMETER ::  nfft =  32  !< maximum length of calls to *fft
@@ -2197,3 +2214,1695 @@
  END SUBROUTINE set99
 
+
+ SUBROUTINE fft991cy_vec( a, work, trigs, ifax, n, isign )
+
+    USE kinds
+
+    IMPLICIT NONE
+
+    REAL(wp),DIMENSION(:,:)     ::  a    !<
+    REAL(wp),DIMENSION(:)       ::  trigs !<
+    REAL(wp),DIMENSION(:,:)     ::  work  !<
+    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) ::  jj     !<
+    INTEGER(iwp) ::  k      !<
+    INTEGER(iwp) ::  la     !<
+    INTEGER(iwp) ::  nb     !<
+    INTEGER(iwp) ::  nblox  !<
+    INTEGER(iwp) ::  nfax   !<
+    INTEGER(iwp) ::  nvex   !<
+    INTEGER(iwp) ::  nx     !<
+    INTEGER(iwp) ::  mm     !<
+
+
+    inc  = 1
+    jump = n
+    lot  = 1
+
+    IF ( ifax(10) /= n )  CALL set99( trigs, ifax, n )
+    nfax = ifax(1)
+    nx   = n + 1
+    IF ( MOD(n,2) == 1 )  nx = n
+    nblox = 1
+    nvex = 1
+
+    IF ( isign == 1 )  THEN
+!
+!--    Backward fft: spectral to gridpoint transform
+
+       istart = 1
+       ia = istart
+       i = istart
+       a(:,i+inc) = 0.5_wp * a(:,i)
+       IF ( MOD(n,2) /= 1 )  THEN
+          i = istart + n * inc
+          a(:,i) = 0.5_wp * a(:,i)
+       ENDIF
+       ia = istart + inc
+       la = 1
+       igo = + 1
+
+       DO  k = 1, nfax
+
+          ifac = ifax(k+1)
+          ierr = -1
+          IF ( igo /= -1 )  THEN
+             CALL rpassm_vec( a(:,ia:), a(:,ia+la*inc:), work, work(:,ifac*la+1:), trigs, inc, 1,  &
+                              n, ifac, la, ierr )
+          ELSE
+             CALL rpassm_vec( work, work(:,la+1:), a(:,ia:), a(:,ia+ifac*la*inc:), trigs, 1, inc,  &
+                              n, ifac, la, ierr )
+          ENDIF
+!
+!--       Following messages shouldn't appear in PALM applications
+          IF ( ierr /= 0 )  THEN
+             SELECT CASE (ierr)
+                CASE (:-1)
+                   WRITE (nout,'(A,I5,A)') ' Vector length =',nvex,', greater than nfft'
+                CASE (0)
+                   WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', not catered for'
+                CASE (1:)
+                   WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', only catered for if la*ifac=n'
+             END SELECT
+             RETURN
+          ENDIF
+
+          la  = ifac * la
+          igo = -igo
+          ia  = istart
+
+       ENDDO
+!
+!--    If necessary, copy results back to a
+       IF ( MOD(nfax,2) /= 0 )  THEN
+          ibase = 1
+          jbase = ia
+          i = ibase
+          j = jbase
+          DO  ii = 1, n
+             a(:,j) = work(:,i)
+             i = i + 1
+             j = j + inc
+          ENDDO
+       ENDIF
+!
+!--    Fill in zeros at end
+       ix = istart + n*inc
+       a(:,ix) = 0.0_wp
+       a(:,ix+inc) = 0.0_wp
+
+    ELSEIF ( isign == -1 )  THEN
+!
+!--    Forward fft: gridpoint to spectral transform
+       istart = 1
+       ia  = istart
+       la  = n
+       igo = + 1
+
+       DO  k = 1, nfax
+
+          ifac = ifax(nfax+2-k)
+          la = la / ifac
+          ierr = -1
+
+          IF ( igo /= -1 )  THEN
+             CALL qpassm_vec( a(:,ia:), a(:,ia+ifac*la*inc:), work, work(:,la+1:), trigs, inc, 1,  &
+                              n, ifac, la, ierr )
+          ELSE
+             CALL qpassm_vec( work, work(:,ifac*la+1:), a(:,ia:), a(:,ia+la*inc:), trigs, 1, inc,  &
+                              n, ifac, la, ierr )
+          ENDIF
+
+!
+!--       Following messages shouldn't appear in PALM applications
+          IF ( ierr /= 0 )  THEN
+             SELECT CASE (ierr)
+                CASE (0)
+                   WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', not catered for'
+                CASE (1:)
+                   WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', only catered for if la*ifac=n'
+             END SELECT
+             RETURN
+          ENDIF
+
+          igo = -igo
+          ia = istart + inc
+
+       ENDDO
+!
+!--    If necessary, copy results back to a
+       IF ( MOD(nfax,2) /= 0 )  THEN
+          ibase = 1
+          jbase = ia
+          i = ibase
+          j = jbase
+          DO  ii = 1, n
+             a(:,j) = work(:,i)
+             i = i + 1
+             j = j + inc
+          ENDDO
+       ENDIF
+!
+!--    Shift a(0) and fill in zero imag parts
+       ix = istart
+       a(:,ix) = a(:,ix+inc)
+       a(:,ix+inc) = 0.0_wp
+
+       IF ( MOD(n,2) /= 1 )  THEN
+          iz = istart + (n+1) * inc
+          a(:,iz) = 0.0_wp
+       ENDIF
+
+    ENDIF
+
+ END SUBROUTINE fft991cy_vec
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> 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)
+!> c       is first real output vector
+!>         equivalence d(1) with c(la*inc2+1)
+!> trigs   is a precalculated list of sines & cosines
+!> inc1    is the addressing increment for a
+!> inc2    is the addressing increment for c
+!> inc3    is the increment between input vectors a
+!> inc4    is the increment between output vectors c
+!> lot     is the number of vectors
+!> n       is the length of the vectors
+!> ifac    is the current factor of n
+!>         la = n/(product of factors used so far)
+!> ierr    is an error indicator:
+!>         0 - pass completed without error
+!>         1 - lot greater than nfft
+!>         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 )
+
+    USE kinds
+
+    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 !<
+
+!
+!-- 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     !<
+
+
+    DATA  sin36/0.587785252292473_wp/, sin72/0.951056516295154_wp/,                                &
+          qrt5/0.559016994374947_wp/,  sin60/0.866025403784437_wp/
+
+
+    ierr = 0
+
+    m = n / ifac
+    iink  = la * inc1
+    jink  = la * inc2
+    ijump = (ifac-1) * iink
+    kstop = ( n-ifac ) / ( 2*ifac )
+
+    ibase = 0
+    jbase = 0
+    igo = ifac - 1
+    IF ( igo == 7 )  igo = 6
+
+    IF (igo < 1 .OR. igo > 6 )  THEN
+       ierr = 2
+       RETURN
+    ENDIF
+
+
+    SELECT CASE ( igo )
+
+!
+!--    Coding for factor 2
+       CASE ( 1 )
+          ia = 1
+          ib = ia + iink
+          ja = 1
+          jb = ja + (2*m-la) * inc2
+
+          IF ( la /= m )  THEN
+
+             DO  l = 1, la
+                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 + ijump
+             ijump = 2 * ijump + iink
+
+             IF ( ja /= jb )  THEN
+
+                DO  k = la, kstop, la
+                   kb = k + k
+                   c1 = trigs(kb+1)
+                   s1 = trigs(kb+2)
+                   jbase = 0
+                   DO  l = 1, la
+                      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)
+                      ibase = ibase + inc1
+                      jbase = jbase + inc2
+                   ENDDO
+
+                   ibase = ibase + ijump
+                   ja = ja + jink
+                   jb = jb - jink
+                ENDDO
+
+                IF ( ja > jb )  RETURN
+
+             ENDIF
+
+             jbase = 0
+             DO  l = 1, la
+                i = ibase
+                j = jbase
+                c(:,ja+j) = a(:,ia+i)
+                d(:,ja+j) = -a(:,ib+i)
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          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))
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          ENDIF
+
+!
+!--    Coding for factor 3
+       CASE ( 2 )
+
+          ia = 1
+          ib = ia + iink
+          ic = ib + iink
+          ja = 1
+          jb = ja + (2*m-la) * inc2
+          jc = jb
+
+          IF ( la /= m )  THEN
+
+             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
+             ibase = ibase + ijump
+             ijump = 2 * ijump + iink
+
+             IF ( ja /= jc )  THEN
+
+                DO  k = la, kstop, la
+                   kb = k + k
+                   kc = kb + kb
+                   c1 = trigs(kb+1)
+                   s1 = trigs(kb+2)
+                   c2 = trigs(kc+1)
+                   s2 = trigs(kc+2)
+
+                   jbase = 0
+                   DO  l = 1, la
+                      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)))
+                         c(ipl,ja+j) = a(ipl,ia+i) + a1
+                         d(ipl,ja+j) = b(ipl,ia+i) + b1
+                         c(ipl,jb+j) = a2 + b3
+                         d(ipl,jb+j) = b2 - a3
+                         c(ipl,jc+j) = a2 - b3
+                         d(ipl,jc+j) = -(b2+a3)
+                      ENDDO
+                      ibase = ibase + inc1
+                      jbase = jbase + inc2
+                   ENDDO
+                   ibase = ibase + ijump
+                   ja = ja + jink
+                   jb = jb + jink
+                   jc = jc - jink
+                ENDDO
+
+                IF ( ja > jc )  RETURN
+
+             ENDIF
+
+             jbase = 0
+             DO  l = 1, la
+                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))
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          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))
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          ENDIF
+
+!
+!--    Coding for factor 4
+       CASE ( 3 )
+
+          ia = 1
+          ib = ia + iink
+          ic = ib + iink
+          id = ic + iink
+          ja = 1
+          jb = ja + (2*m-la) * inc2
+          jc = jb + 2*m*inc2
+          jd = jb
+
+          IF ( la /= m )  THEN
+
+             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))
+                c(:,jb+j) = a(:,ia+i) - a(:,ic+i)
+                d(:,jb+j) = a(:,id+i) - a(:,ib+i)
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+             ja = ja + jink
+             jink = 2 * jink
+             jb = jb + jink
+             jc = jc - jink
+             jd = jd - jink
+             ibase = ibase + ijump
+             ijump = 2 * ijump + iink
+
+             IF ( jb /= jc )  THEN
+
+                DO  k = la, kstop, la
+                   kb = k + k
+                   kc = kb + kb
+                   kd = kc + kb
+                   c1 = trigs(kb+1)
+                   s1 = trigs(kb+2)
+                   c2 = trigs(kc+1)
+                   s2 = trigs(kc+2)
+                   c3 = trigs(kd+1)
+                   s3 = trigs(kd+2)
+                   jbase = 0
+                   DO  l = 1, la
+                      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))
+                         c(ipl,ja+j) = a0 + a1
+                         c(ipl,jc+j) = a0 - a1
+                         d(ipl,ja+j) = b0 + b1
+                         d(ipl,jc+j) = b1 - b0
+                         c(ipl,jb+j) = a2 + b3
+                         c(ipl,jd+j) = a2 - b3
+                         d(ipl,jb+j) = b2 - a3
+                         d(ipl,jd+j) = -(b2+a3)
+                      ENDDO
+                      ibase = ibase + inc1
+                      jbase = jbase + inc2
+                   ENDDO
+                   ibase = ibase + ijump
+                   ja = ja + jink
+                   jb = jb + jink
+                   jc = jc - jink
+                   jd = jd - jink
+                ENDDO
+
+                IF ( jb > jc )  RETURN
+
+             ENDIF
+
+             sin45 = SQRT( 0.5_wp )
+             jbase = 0
+             DO  l = 1, la
+                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))
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          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))
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          ENDIF
+
+!
+!--    Coding for factor 5
+       CASE ( 4 )
+
+          ia = 1
+          ib = ia + iink
+          ic = ib + iink
+          id = ic + iink
+          ie = id + iink
+          ja = 1
+          jb = ja + (2*m-la) * inc2
+          jc = jb + 2*m*inc2
+          jd = jc
+          je = jb
+
+          IF ( la /= m )  THEN
+
+             DO  l = 1, la
+                i = ibase
+                j = jbase
+                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)
+                   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
+             ibase = ibase + ijump
+             ijump = 2 * ijump + iink
+
+             IF ( jb /= jd )  THEN
+
+                DO  k = la, kstop, la
+                   kb = k + k
+                   kc = kb + kb
+                   kd = kc + kb
+                   ke = kd + kb
+                   c1 = trigs(kb+1)
+                   s1 = trigs(kb+2)
+                   c2 = trigs(kc+1)
+                   s2 = trigs(kc+2)
+                   c3 = trigs(kd+1)
+                   s3 = trigs(kd+2)
+                   c4 = trigs(ke+1)
+                   s4 = trigs(ke+2)
+                   jbase = 0
+                   DO  l = 1, la
+                      i = ibase
+                      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)
+                         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)
+                         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,jb+j) = b10 - b11
+                         d(ipl,je+j) = -(b10+b11)
+                         d(ipl,jc+j) = b20 - b21
+                         d(ipl,jd+j) = -(b20+b21)
+                      ENDDO
+                      ibase = ibase + inc1
+                      jbase = jbase + inc2
+                   ENDDO
+                   ibase = ibase + ijump
+                   ja = ja + jink
+                   jb = jb + jink
+                   jc = jc + jink
+                   jd = jd - jink
+                   je = je - jink
+                ENDDO
+
+                IF ( jb > jd )  RETURN
+
+             ENDIF
+
+             jbase = 0
+             DO  l = 1, la
+                i = ibase
+                j = jbase
+                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)
+                   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
+                ENDDO
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          ELSE
+
+             z = 1.0_wp / REAL( n, KIND=wp )
+             zqrt5  = z * qrt5
+             zsin36 = z * sin36
+             zsin72 = z * sin72
+             DO  l = 1, la
+                i = ibase
+                j = jbase
+                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))
+                   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
+                ENDDO
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          ENDIF
+
+!
+!--    Coding for factor 6
+       CASE ( 5 )
+
+          ia = 1
+          ib = ia + iink
+          ic = ib + iink
+          id = ic + iink
+          ie = id + iink
+          if = ie + iink
+          ja = 1
+          jb = ja + (2*m-la) * inc2
+          jc = jb + 2*m*inc2
+          jd = jc + 2*m*inc2
+          je = jc
+          jf = jb
+
+          IF ( la /= m )  THEN
+
+             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) = (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
+                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
+             ibase = ibase + ijump
+             ijump = 2 * ijump + iink
+
+             IF ( jc /= jd )  THEN
+
+                DO  k = la, kstop, la
+                   kb = k + k
+                   kc = kb + kb
+                   kd = kc + kb
+                   ke = kd + kb
+                   kf = ke + kb
+                   c1 = trigs(kb+1)
+                   s1 = trigs(kb+2)
+                   c2 = trigs(kc+1)
+                   s2 = trigs(kc+2)
+                   c3 = trigs(kd+1)
+                   s3 = trigs(kd+2)
+                   c4 = trigs(ke+1)
+                   s4 = trigs(ke+2)
+                   c5 = trigs(kf+1)
+                   s5 = trigs(kf+2)
+                   jbase = 0
+                   DO  l = 1, la
+                      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
+                         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))
+                         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,jf+j) = a20 + b21
+                         d(ipl,jf+j) = a21 + b20
+                      ENDDO
+                      ibase = ibase + inc1
+                      jbase = jbase + inc2
+                   ENDDO
+                   ibase = ibase + ijump
+                   ja = ja + jink
+                   jb = jb + jink
+                   jc = jc + jink
+                   jd = jd - jink
+                   je = je - jink
+                   jf = jf - jink
+                ENDDO
+
+                IF ( jc > jd )  RETURN
+
+             ENDIF
+
+             jbase = 0
+             DO  l = 1, la
+                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))
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          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)
+                ENDDO
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          ENDIF
+
+!
+!--    Coding for factor 8
+       CASE ( 6 )
+
+          IF ( la /= m )  THEN
+             ierr = 3
+             RETURN
+          ENDIF
+
+          ia = 1
+          ib = ia + iink
+          ic = ib + iink
+          id = ic + iink
+          ie = id + iink
+          if = ie + iink
+          ig = if + iink
+          ih = ig + iink
+          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 )
+          zsin45 = z * SQRT( 0.5_wp )
+
+          DO  l = 1, la
+             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))
+             ibase = ibase + inc1
+             jbase = jbase + inc2
+          ENDDO
+
+    END SELECT
+
+ 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 )
+
+    USE kinds
+
+    IMPLICIT NONE
+
+    INTEGER(iwp) ::  ierr !<
+    INTEGER(iwp) ::  ifac !<
+    INTEGER(iwp) ::  inc1 !<
+    INTEGER(iwp) ::  inc2 !<
+    INTEGER(iwp) ::  la   !<
+    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       !<
+
+
+    DATA  sin36/0.587785252292473_wp/, sin72/0.951056516295154_wp/,                                &
+          qrt5/0.559016994374947_wp/,  sin60/0.866025403784437_wp/
+
+
+    ierr = 0
+
+    m = n / ifac
+    iink = la * inc1
+    jink = la * inc2
+    jump = (ifac-1) * jink
+    kstop = (n-ifac) / (2*ifac)
+
+    ibase = 0
+    jbase = 0
+    igo = ifac - 1
+    IF ( igo == 7 )  igo = 6
+    IF ( igo < 1  .OR.  igo > 6 )  THEN
+       ierr = 2
+       RETURN
+    ENDIF
+
+
+    SELECT CASE ( igo )
+!
+!--    Coding for factor 2
+       CASE ( 1 )
+
+          ia = 1
+          ib = ia + (2*m-la) * inc1
+          ja = 1
+          jb = ja + jink
+
+          IF ( la /= m )  THEN
+
+             DO  l = 1, la
+                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
+
+             ia = ia + iink
+             iink = 2 * iink
+             ib = ib - iink
+             ibase = 0
+             jbase = jbase + jump
+             jump = 2 * jump + jink
+
+             IF ( ia /= ib )  THEN
+
+                DO  k = la, kstop, la
+                   kb = k + k
+                   c1 = trigs(kb+1)
+                   s1 = trigs(kb+2)
+                   ibase = 0
+                   DO  l = 1, la
+                      i = ibase
+                      j = jbase
+                      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))
+                      ibase = ibase + inc1
+                      jbase = jbase + inc2
+                   ENDDO
+                   ia = ia + iink
+                   ib = ib - iink
+                   jbase = jbase + jump
+                ENDDO
+
+                IF ( ia > ib )  RETURN
+
+             ENDIF
+
+             ibase = 0
+             DO  l = 1, la
+                i = ibase
+                j = jbase
+                c(:,ja+j) = a(:,ia+i)
+                c(:,jb+j) = -b(:,ia+i)
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          ELSE
+
+             DO  l = 1, la
+                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))
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          ENDIF
+
+!
+!--    Coding for factor 3
+       CASE ( 2 )
+
+          ia = 1
+          ib = ia + (2*m-la) * inc1
+          ic = ib
+          ja = 1
+          jb = ja + jink
+          jc = jb + jink
+
+          IF ( la /= m )  THEN
+
+             DO  l = 1, la
+                i = ibase
+                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)))
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+             ia = ia + iink
+             iink = 2 * iink
+             ib = ib + iink
+             ic = ic - iink
+             jbase = jbase + jump
+             jump = 2 * jump + jink
+
+             IF ( ia /= ic )  THEN
+
+                DO  k = la, kstop, la
+                   kb = k + k
+                   kc = kb + kb
+                   c1 = trigs(kb+1)
+                   s1 = trigs(kb+2)
+                   c2 = trigs(kc+1)
+                   s2 = trigs(kc+2)
+                   ibase = 0
+                   DO  l = 1, la
+                      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))))
+                      ibase = ibase + inc1
+                      jbase = jbase + inc2
+                   ENDDO
+                   ia = ia + iink
+                   ib = ib + iink
+                   ic = ic - iink
+                   jbase = jbase + jump
+                ENDDO
+
+                IF ( ia > ic )  RETURN
+
+             ENDIF
+
+             ibase = 0
+             DO  l = 1, la
+                i = ibase
+                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))
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          ELSE
+
+             ssin60 = 2.0_wp * sin60
+             DO  l = 1, la
+                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))
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          ENDIF
+
+!
+!--    Coding for factor 4
+       CASE ( 3 )
+
+          ia = 1
+          ib = ia + (2*m-la) * inc1
+          ic = ib + 2*m*inc1
+          id = ib
+          ja = 1
+          jb = ja + jink
+          jc = jb + jink
+          jd = jc + jink
+
+          IF ( la /= m )  THEN
+
+             DO  l = 1, la
+                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)
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+             ia = ia + iink
+             iink = 2 * iink
+             ib = ib + iink
+             ic = ic - iink
+             id = id - iink
+             jbase = jbase + jump
+             jump = 2 * jump + jink
+
+             IF ( ib /= ic )  THEN
+
+                DO  k = la, kstop, la
+                   kb = k + k
+                   kc = kb + kb
+                   kd = kc + kb
+                   c1 = trigs(kb+1)
+                   s1 = trigs(kb+2)
+                   c2 = trigs(kc+1)
+                   s2 = trigs(kc+2)
+                   c3 = trigs(kd+1)
+                   s3 = trigs(kd+2)
+                   ibase = 0
+                   DO  l = 1, la
+                      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)))
+                      ibase = ibase + inc1
+                      jbase = jbase + inc2
+                   ENDDO
+                   ia = ia + iink
+                   ib = ib + iink
+                   ic = ic - iink
+                   id = id - iink
+                   jbase = jbase + jump
+                ENDDO
+
+                IF ( ib > ic )  RETURN
+
+             ENDIF
+
+             ibase = 0
+             sin45 = SQRT( 0.5_wp )
+             DO  l = 1, la
+                i = ibase
+                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(:,jc+j) = b(:,ib+i) - b(:,ia+i)
+                c(:,jd+j) = -sin45*((a(:,ia+i)-a(:,ib+i))+(b(:,ia+i)+b(:,ib+i)))
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          ELSE
+
+             DO  l = 1, la
+                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))
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          ENDIF
+
+!
+!--    Coding for factor 5
+       CASE ( 4 )
+
+          ia = 1
+          ib = ia + (2*m-la) * inc1
+          ic = ib + 2*m*inc1
+          id = ic
+          ie = ib
+          ja = 1
+          jb = ja + jink
+          jc = jb + jink
+          jd = jc + jink
+          je = jd + jink
+
+          IF ( la /= m )  THEN
+
+             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.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
+             ENDDO
+             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
+
+             IF ( ib /= id )  THEN
+
+                DO  k = la, kstop, la
+                   kb = k + k
+                   kc = kb + kb
+                   kd = kc + kb
+                   ke = kd + kb
+                   c1 = trigs(kb+1)
+                   s1 = trigs(kb+2)
+                   c2 = trigs(kc+1)
+                   s2 = trigs(kc+2)
+                   c3 = trigs(kd+1)
+                   s3 = trigs(kd+2)
+                   c4 = trigs(ke+1)
+                   s4 = trigs(ke+2)
+                   ibase = 0
+                   DO  l = 1, la
+                      i = ibase
+                      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     )
+                      ENDDO
+                      ibase = ibase + inc1
+                      jbase = jbase + inc2
+                   ENDDO
+                   ia = ia + iink
+                   ib = ib + iink
+                   ic = ic + iink
+                   id = id - iink
+                   ie = ie - iink
+                   jbase = jbase + jump
+                ENDDO
+
+                IF ( ib > id )  RETURN
+
+             ENDIF
+
+             ibase = 0
+             DO  l = 1, la
+                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))
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          ELSE
+
+             qqrt5  = 2.0_wp * qrt5
+             ssin36 = 2.0_wp * sin36
+             ssin72 = 2.0_wp * sin72
+             DO  l = 1, la
+                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))
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          ENDIF
+
+!
+!--    Coding for factor 6
+       CASE ( 5 )
+
+          ia = 1
+          ib = ia + (2*m-la) * inc1
+          ic = ib + 2*m*inc1
+          id = ic + 2*m*inc1
+          ie = ic
+          if = ib
+          ja = 1
+          jb = ja + jink
+          jc = jb + jink
+          jd = jc + jink
+          je = jd + jink
+          jf = je + jink
+
+          IF ( la /= m )  THEN
+
+             DO  l = 1, la
+                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)))
+                ibase = ibase + inc1
+                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
+             jbase = jbase + jump
+             jump = 2 * jump + jink
+
+             IF ( ic /= id )  THEN
+
+                DO  k = la, kstop, la
+                   kb = k + k
+                   kc = kb + kb
+                   kd = kc + kb
+                   ke = kd + kb
+                   kf = ke + kb
+                   c1 = trigs(kb+1)
+                   s1 = trigs(kb+2)
+                   c2 = trigs(kc+1)
+                   s2 = trigs(kc+2)
+                   c3 = trigs(kd+1)
+                   s3 = trigs(kd+2)
+                   c4 = trigs(ke+1)
+                   s4 = trigs(ke+2)
+                   c5 = trigs(kf+1)
+                   s5 = trigs(kf+2)
+                   ibase = 0
+                   DO  l = 1, la
+                      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     )
+
+                      ENDDO
+                      ibase = ibase + inc1
+                      jbase = jbase + inc2
+                   ENDDO
+                   ia = ia + iink
+                   ib = ib + iink
+                   ic = ic + iink
+                   id = id - iink
+                   ie = ie - iink
+                   if = if - iink
+                   jbase = jbase + jump
+                ENDDO
+
+                IF ( ic > id )  RETURN
+
+             ENDIF
+
+             ibase = 0
+             DO  l = 1, la
+                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))
+                ENDDO
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          ELSE
+
+             ssin60 = 2.0_wp * sin60
+             DO  l = 1, la
+                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)))
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          ENDIF
+
+!
+!--    Coding for factor 8
+       CASE ( 6 )
+
+          IF ( la /= m )  THEN
+             ierr = 3
+             RETURN
+          ENDIF
+
+          ia = 1
+          ib = ia + la*inc1
+          ic = ib + 2*la*inc1
+          id = ic + 2*la*inc1
+          ie = id + 2*la*inc1
+          ja = 1
+          jb = ja + jink
+          jc = jb + jink
+          jd = jc + jink
+          je = jd + jink
+          jf = je + jink
+          jg = jf + jink
+          jh = jg + jink
+          ssin45 = SQRT( 2.0_wp )
+
+          DO  l = 1, la
+             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)))
+             ibase = ibase + inc1
+             jbase = jbase + inc2
+          ENDDO
+
+    END SELECT
+
+ END SUBROUTINE rpassm_vec
+
  END MODULE temperton_fft
Index: /palm/trunk/SOURCE/transpose.f90
===================================================================
--- /palm/trunk/SOURCE/transpose.f90	(revision 4365)
+++ /palm/trunk/SOURCE/transpose.f90	(revision 4366)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! modifications for NEC vectorization
+! 
+! 4360 2020-01-07 11:25:50Z suehring
 ! Added missing OpenMP directives
 ! 
@@ -266,4 +269,7 @@
         ONLY:  cpu_log, cpu_log_nowait, log_point_s
 
+    USE fft_xy,                                                                &
+        ONLY:  f_vec, temperton_fft_vec
+
     USE indices,                                                               &
         ONLY:  nnx, nx, nxl, nxr, nyn, nys, nz
@@ -282,4 +288,5 @@
     INTEGER(iwp) ::  k  !<
     INTEGER(iwp) ::  l  !<
+    INTEGER(iwp) ::  mm !<
     INTEGER(iwp) ::  xs !<
 
@@ -292,31 +299,49 @@
 #endif
 
-
-!
-!-- If the PE grid is one-dimensional along y, the array has only to be
-!-- reordered locally and therefore no transposition has to be done.
+    !
+    !-- If the PE grid is one-dimensional along y, the array has only to be
+    !-- reordered locally and therefore no transposition has to be done.
     IF ( pdims(1) /= 1 )  THEN
 
 #if defined( __parallel )
 !
-!--    Reorder input array for transposition
-!$OMP  PARALLEL PRIVATE ( i, j, k, l, xs )
-       DO  l = 0, pdims(1) - 1
-          xs = 0 + l * nnx
-#if __acc_fft_device
-          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
-          !$ACC PRESENT(work, f_in)
-#endif
-          !$OMP DO
-          DO  k = nzb_x, nzt_x
-             DO  i = xs, xs + nnx - 1
-                DO  j = nys_x, nyn_x
-                   work(j,i-xs+1,k,l) = f_in(i,j,k)
+!--    Reorder input array for transposition. Data from the vectorized Temperton-fft is stored in
+!--    different array format (f_vec).
+       IF ( temperton_fft_vec )  THEN
+
+          DO  l = 0, pdims(1) - 1
+             xs = 0 + l * nnx
+             DO  k = nzb_x, nzt_x
+                DO  i = xs, xs + nnx - 1
+                   DO  j = nys_x, nyn_x
+                      mm = j-nys_x+1+(k-nzb_x)*(nyn_x-nys_x+1)
+                      work(j,i-xs+1,k,l) = f_vec(mm,i)
+                   ENDDO
                 ENDDO
              ENDDO
           ENDDO
-          !$OMP END DO NOWAIT
-       ENDDO
-!$OMP  END PARALLEL
+
+       ELSE
+
+          !$OMP  PARALLEL PRIVATE ( i, j, k, l, xs )
+          DO  l = 0, pdims(1) - 1
+             xs = 0 + l * nnx
+#if __acc_fft_device
+             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+             !$ACC PRESENT(work, f_in)
+#endif
+             !$OMP DO
+             DO  k = nzb_x, nzt_x
+                DO  i = xs, xs + nnx - 1
+                   DO  j = nys_x, nyn_x
+                      work(j,i-xs+1,k,l) = f_in(i,j,k)
+                   ENDDO
+                ENDDO
+             ENDDO
+             !$OMP END DO NOWAIT
+          ENDDO
+          !$OMP  END PARALLEL
+
+       ENDIF
 
 !
@@ -836,4 +861,7 @@
         ONLY:  cpu_log, cpu_log_nowait, log_point_s
 
+    USE fft_xy,                                                                &
+        ONLY:  f_vec, temperton_fft_vec
+
     USE indices,                                                               &
         ONLY:  nnx, nx, nxl, nxr, nyn, nys, nz
@@ -852,4 +880,5 @@
     INTEGER(iwp) ::  k  !<
     INTEGER(iwp) ::  l  !<
+    INTEGER(iwp) ::  mm !<
     INTEGER(iwp) ::  xs !<
 
@@ -914,24 +943,44 @@
 
 !
-!--    Reorder transposed array
-!$OMP  PARALLEL PRIVATE ( i, j, k, l, xs )
-       DO  l = 0, pdims(1) - 1
-          xs = 0 + l * nnx
-#if __acc_fft_device
-          !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
-          !$ACC PRESENT(f_out, work)
-#endif
-          !$OMP DO
-          DO  k = nzb_x, nzt_x
-             DO  i = xs, xs + nnx - 1
-                DO  j = nys_x, nyn_x
-                   f_out(i,j,k) = work(j,i-xs+1,k,l)
+!--    Reorder transposed array.
+!--    Data for the vectorized Temperton-fft is stored in different array format (f_vec) which saves
+!--    additional data copy in fft_x.
+       IF ( temperton_fft_vec )  THEN
+
+          DO  l = 0, pdims(1) - 1
+             xs = 0 + l * nnx
+             DO  k = nzb_x, nzt_x
+                DO  i = xs, xs + nnx - 1
+                   DO  j = nys_x, nyn_x
+                      mm = j-nys_x+1+(k-nzb_x)*(nyn_x-nys_x+1)
+                      f_vec(mm,i) = work(j,i-xs+1,k,l)
+                   ENDDO
                 ENDDO
              ENDDO
           ENDDO
-          !$OMP END DO NOWAIT
-       ENDDO
-!$OMP  END PARALLEL
-#endif
+
+       ELSE
+
+          !$OMP  PARALLEL PRIVATE ( i, j, k, l, xs )
+          DO  l = 0, pdims(1) - 1
+             xs = 0 + l * nnx
+#if __acc_fft_device
+             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i,j,k) &
+             !$ACC PRESENT(f_out, work)
+#endif
+             !$OMP DO
+             DO  k = nzb_x, nzt_x
+                DO  i = xs, xs + nnx - 1
+                   DO  j = nys_x, nyn_x
+                      f_out(i,j,k) = work(j,i-xs+1,k,l)
+                   ENDDO
+                ENDDO
+             ENDDO
+             !$OMP END DO NOWAIT
+          ENDDO
+          !$OMP  END PARALLEL
+#endif
+
+       ENDIF
 
     ENDIF
