Index: palm/trunk/SCRIPTS/palmbuild
===================================================================
--- palm/trunk/SCRIPTS/palmbuild	(revision 3724)
+++ palm/trunk/SCRIPTS/palmbuild	(revision 3725)
@@ -27,4 +27,8 @@
 # -----------------
 # $Id$
+# ssh-calls for compilations on remote systems modified to avoid output
+# of login messages on specific systems changed again (reverted as before r3549)
+# 
+# 3549 2018-11-21 15:44:44Z raasch
 # ssh-calls for compilations on remote systems modified to avoid output
 # of login messages on specific systems
@@ -571,6 +575,6 @@
        fi
        make_call_string="make  -f Makefile_utilities  $make_options  F90=$compiler_name  F90_SER=$compiler_name_ser  COPT=\"$cpp_options\"  F90FLAGS=\"$compiler_options\"  LDFLAGS=\"$linker_options\" "
-       echo "$login_init_cmd $module_commands cd ${make_depository}; echo $make_call_string > LAST_MAKE_CALL; chmod u+x LAST_MAKE_CALL; $make_call_string; [[ \$? != 0 ]] && echo MAKE_ERROR" | ssh  -q  $ssh_key  ${remote_username}@${remote_ip} 2>/dev/null | tee ${configuration_identifier}_last_make_protocol
-###       ssh  -q  $ssh_key  ${remote_username}@${remote_ip}  "$login_init_cmd $module_commands cd ${make_depository}; echo $make_call_string > LAST_MAKE_CALL; chmod u+x LAST_MAKE_CALL; $make_call_string; [[ \$? != 0 ]] && echo MAKE_ERROR"  2>&1 | tee ${configuration_identifier}_last_make_protocol
+###       echo "$login_init_cmd $module_commands cd ${make_depository}; echo $make_call_string > LAST_MAKE_CALL; chmod u+x LAST_MAKE_CALL; $make_call_string; [[ \$? != 0 ]] && echo MAKE_ERROR" | ssh  -q  $ssh_key  ${remote_username}@${remote_ip} 2>/dev/null | tee ${configuration_identifier}_last_make_protocol
+       ssh  -q  $ssh_key  ${remote_username}@${remote_ip}  "$login_init_cmd $module_commands cd ${make_depository}; echo $make_call_string > LAST_MAKE_CALL; chmod u+x LAST_MAKE_CALL; $make_call_string; [[ \$? != 0 ]] && echo MAKE_ERROR"  2>&1 | tee ${configuration_identifier}_last_make_protocol
 
        if [[ $(grep -c MAKE_ERROR ${configuration_identifier}_last_make_protocol) != 0 ]]
@@ -611,6 +615,6 @@
        fi
        make_call_string="make  $make_options  PROG=$program_name  F90=$compiler_name  COPT=\"$cpp_options\"  F90FLAGS=\"$compiler_options\"  LDFLAGS=\"$linker_options\" "
-       echo "$login_init_cmd $module_commands cd ${make_depository}; echo $make_call_string > LAST_MAKE_CALL; chmod u+x LAST_MAKE_CALL; $make_call_string; [[ \$? != 0 ]] && echo MAKE_ERROR" | ssh  -q  $ssh_key  ${remote_username}@${remote_ip} 2>&1 | tee ${configuration_identifier}_last_make_protocol
-###       ssh  -q  $ssh_key  ${remote_username}@${remote_ip}  "$login_init_cmd $module_commands cd ${make_depository}; echo $make_call_string > LAST_MAKE_CALL; chmod u+x LAST_MAKE_CALL; $make_call_string; [[ \$? != 0 ]] && echo MAKE_ERROR"  2>&1 | tee ${configuration_identifier}_last_make_protocol
+###       echo "$login_init_cmd $module_commands cd ${make_depository}; echo $make_call_string > LAST_MAKE_CALL; chmod u+x LAST_MAKE_CALL; $make_call_string; [[ \$? != 0 ]] && echo MAKE_ERROR" | ssh  -q  $ssh_key  ${remote_username}@${remote_ip} 2>&1 | tee ${configuration_identifier}_last_make_protocol
+       ssh  -q  $ssh_key  ${remote_username}@${remote_ip}  "$login_init_cmd $module_commands cd ${make_depository}; echo $make_call_string > LAST_MAKE_CALL; chmod u+x LAST_MAKE_CALL; $make_call_string; [[ \$? != 0 ]] && echo MAKE_ERROR"  2>&1 | tee ${configuration_identifier}_last_make_protocol
 
        if [[ $(grep -c MAKE_ERROR ${configuration_identifier}_last_make_protocol) != 0 ]]
Index: palm/trunk/SCRIPTS/palmrun
===================================================================
--- palm/trunk/SCRIPTS/palmrun	(revision 3724)
+++ palm/trunk/SCRIPTS/palmrun	(revision 3725)
@@ -27,4 +27,7 @@
 # -----------------
 # $Id$
+# error messages for failed restarts extended
+# 
+# 3682 2019-01-18 17:01:54Z knoop
 # ssh-call for submitting batch jobs on remote systems modified again to avoid output
 # of login messages on specific systems
@@ -2799,5 +2802,5 @@
 
           # CHECK, IF RESTART JOB HAS BEEN STARTED
-       if [[ $( grep -c "+++ palmrun killed"  palmrun_restart.log ) != 0 ]]
+       if [[ $( grep -c "+++ palmrun crashed"  palmrun_restart.log ) != 0 ]]
        then
           printf "\n$dashes\n  +++ creating restart run failed \n"
@@ -2805,4 +2808,10 @@
           rm  palmrun_restart.log
           exit
+       elif [[ $( grep -c "*** palmrun finished"  palmrun_restart.log ) != 1 ]]
+       then
+          printf "\n$dashes\n  +++ creating restart run failed, probably due to network problems\n"
+          locat=create_restart
+          rm  palmrun_restart.log
+          exit       
        else
           printf "\n$dashes\n  *** restart run initiated \n"
Index: palm/trunk/SOURCE/buoyancy.f90
===================================================================
--- palm/trunk/SOURCE/buoyancy.f90	(revision 3724)
+++ palm/trunk/SOURCE/buoyancy.f90	(revision 3725)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! unused variables removed
+! 
+! 3655 2019-01-07 16:51:22Z knoop
 ! nopointer option removed
 ! 
@@ -160,6 +163,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nzb,       &
-                  nzt
+           ONLY:  nxl, nxlu, nxr, nyn, nys, nzb, nzt
 
 
@@ -251,5 +253,5 @@
 
        USE indices,                                                            &
-           ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt
+           ONLY:  nzb, nzt
 
        IMPLICIT NONE
Index: palm/trunk/SOURCE/data_log.f90
===================================================================
--- palm/trunk/SOURCE/data_log.f90	(revision 3724)
+++ palm/trunk/SOURCE/data_log.f90	(revision 3725)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! preprocessor directives removed to avoid compiler warnings
+! 
+! 3655 2019-01-07 16:51:22Z knoop
 ! Corrected "Former revisions" section
 ! 
@@ -61,6 +64,4 @@
  SUBROUTINE data_log( array, i1, i2, j1, j2, k1, k2 )
  
-#if defined( __logging )
-
     USE control_parameters,                                                    &
         ONLY:  log_message, simulated_time
@@ -98,5 +99,4 @@
     WRITE ( 20 )  array
 
-#endif
  END SUBROUTINE data_log
 
@@ -110,6 +110,4 @@
  
  SUBROUTINE data_log_2d( array, i1, i2, j1, j2)
-
-#if defined( __logging )
 
     USE control_parameters,                                                    &
@@ -146,5 +144,4 @@
     WRITE ( 20 )  array
 
-#endif
  END SUBROUTINE data_log_2d
 
@@ -158,6 +155,4 @@
  
  SUBROUTINE data_log_2d_int( array, i1, i2, j1, j2)
-
-#if defined( __logging )
 
     USE control_parameters,                                                    &
@@ -194,4 +189,3 @@
     WRITE ( 20 )  array
 
-#endif
  END SUBROUTINE data_log_2d_int
Index: palm/trunk/SOURCE/gust_mod.f90
===================================================================
--- palm/trunk/SOURCE/gust_mod.f90	(revision 3724)
+++ palm/trunk/SOURCE/gust_mod.f90	(revision 3725)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! dummy statement modified to avoid compiler warnings about unused variables
+! 
+! 3685 2019-01-21 01:02:11Z knoop
 ! Some interface calls moved to module_interface + cleanup
 ! 
@@ -445,5 +448,7 @@
 !
 !--    Next line is just to avoid compiler warnings about unused variables. You may remove it.
-       IF ( found .AND. two_d )  idum = av + LEN( variable ) + LEN( grid ) + local_pf(nxl,nys,nzb_do) + fill_value
+       IF ( found .AND. two_d )  THEN
+          idum = av + LEN( variable ) + LEN( grid // mode ) + local_pf(nxl,nys,nzb_do) + fill_value
+       ENDIF
 
     END SUBROUTINE gust_data_output_2d
Index: palm/trunk/SOURCE/poismg_mod.f90
===================================================================
--- palm/trunk/SOURCE/poismg_mod.f90	(revision 3724)
+++ palm/trunk/SOURCE/poismg_mod.f90	(revision 3725)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! unused subroutine removed
+! 
+! 3655 2019-01-07 16:51:22Z knoop
 ! unnecessary check eliminated
 ! 
@@ -143,5 +146,5 @@
     INTERFACE sort_k_to_even_odd_blocks
        MODULE PROCEDURE sort_k_to_even_odd_blocks
-       MODULE PROCEDURE sort_k_to_even_odd_blocks_int
+!       MODULE PROCEDURE sort_k_to_even_odd_blocks_int
        MODULE PROCEDURE sort_k_to_even_odd_blocks_1d
     END INTERFACE sort_k_to_even_odd_blocks
@@ -1152,58 +1155,58 @@
 !> Version for 2D-INTEGER arrays
 !------------------------------------------------------------------------------!
-    SUBROUTINE sort_k_to_even_odd_blocks_int( i_mg , glevel )
-
-
-       USE indices,                                                            &
-           ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
-
-       IMPLICIT NONE
-
-       INTEGER(iwp), INTENT(IN) ::  glevel  !< grid level
-
-       INTEGER(iwp), DIMENSION(nzb:nzt_mg(glevel)+1,                           &
-                               nys_mg(glevel)-1:nyn_mg(glevel)+1,              &
-                               nxl_mg(glevel)-1:nxr_mg(glevel)+1) ::           &
-                                    i_mg    !< array to be sorted 
+!    SUBROUTINE sort_k_to_even_odd_blocks_int( i_mg , glevel )
+!
+!
+!       USE indices,                                                            &
+!           ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
+!
+!       IMPLICIT NONE
+!
+!       INTEGER(iwp), INTENT(IN) ::  glevel  !< grid level
+!
+!       INTEGER(iwp), DIMENSION(nzb:nzt_mg(glevel)+1,                           &
+!                               nys_mg(glevel)-1:nyn_mg(glevel)+1,              &
+!                               nxl_mg(glevel)-1:nxr_mg(glevel)+1) ::           &
+!                                    i_mg    !< array to be sorted
 !
 !--    Local variables
-       INTEGER(iwp) :: i        !< index variabel along x
-       INTEGER(iwp) :: j        !< index variable along y
-       INTEGER(iwp) :: k        !< index variable along z
-       INTEGER(iwp) :: l        !< grid level 
-       INTEGER(iwp) :: ind      !< index variable along z
-       INTEGER(iwp),DIMENSION(nzb:nzt_mg(glevel)+1) ::  tmp  !< temporary odd-even sorted array
-
-
-       CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'start' )
-
-       l = glevel
-       ind_even_odd = even_odd_level(l)
-
-       DO  i = nxl_mg(l)-1, nxr_mg(l)+1
-          DO  j = nys_mg(l)-1, nyn_mg(l)+1
-
+!       INTEGER(iwp) :: i        !< index variabel along x
+!       INTEGER(iwp) :: j        !< index variable along y
+!       INTEGER(iwp) :: k        !< index variable along z
+!       INTEGER(iwp) :: l        !< grid level
+!       INTEGER(iwp) :: ind      !< index variable along z
+!       INTEGER(iwp),DIMENSION(nzb:nzt_mg(glevel)+1) ::  tmp  !< temporary odd-even sorted array
+!
+!
+!       CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'start' )
+!
+!       l = glevel
+!       ind_even_odd = even_odd_level(l)
+!
+!       DO  i = nxl_mg(l)-1, nxr_mg(l)+1
+!          DO  j = nys_mg(l)-1, nyn_mg(l)+1
+!
 !
 !--          Sort the data with even k index
-             ind = nzb-1
-             DO  k = nzb, nzt_mg(l), 2
-                ind = ind + 1
-                tmp(ind) = i_mg(k,j,i)
-             ENDDO
+!             ind = nzb-1
+!             DO  k = nzb, nzt_mg(l), 2
+!                ind = ind + 1
+!                tmp(ind) = i_mg(k,j,i)
+!             ENDDO
 !
 !--          Sort the data with odd k index
-             DO  k = nzb+1, nzt_mg(l)+1, 2
-                ind = ind + 1
-                tmp(ind) = i_mg(k,j,i)
-             ENDDO
-
-             i_mg(:,j,i) = tmp
-
-          ENDDO
-       ENDDO
-
-       CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'stop' )
-
-    END SUBROUTINE sort_k_to_even_odd_blocks_int
+!             DO  k = nzb+1, nzt_mg(l)+1, 2
+!                ind = ind + 1
+!                tmp(ind) = i_mg(k,j,i)
+!             ENDDO
+!
+!             i_mg(:,j,i) = tmp
+!
+!          ENDDO
+!       ENDDO
+!
+!       CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'stop' )
+!
+!    END SUBROUTINE sort_k_to_even_odd_blocks_int
 
 
Index: palm/trunk/SOURCE/temperton_fft_mod.f90
===================================================================
--- palm/trunk/SOURCE/temperton_fft_mod.f90	(revision 3724)
+++ palm/trunk/SOURCE/temperton_fft_mod.f90	(revision 3725)
@@ -9,4 +9,7 @@
 ! -----------------
 ! $Id$
+! GOTO statements replaced, file re-formatted corresponding to coding standards
+! 
+! 3241 2018-09-12 15:02:00Z raasch
 ! unused variables removed
 ! 
@@ -25,5 +28,4 @@
 ! Module renamed
 ! 
-! 
 ! 1682 2015-10-07 23:56:08Z knoop
 ! Code annotations made doxygen readable 
@@ -70,5 +72,5 @@
 ! Description:
 ! ------------
-!> Calls fortran-versions of fft's.
+!> Calls Fortran-versions of fft's.
 !> 
 !> Method:
@@ -115,5 +117,5 @@
 !> dimension a(n),work(n),trigs(n),ifax(1)
 !------------------------------------------------------------------------------!
-  SUBROUTINE fft991cy(a,work,trigs,ifax,inc,jump,n,lot,isign)
+ SUBROUTINE fft991cy( a, work, trigs, ifax, inc, jump, n, lot, isign )
 
     USE kinds
@@ -121,5 +123,4 @@
     IMPLICIT NONE
 
-    !  Scalar arguments 
     INTEGER(iwp) ::  inc   !< 
     INTEGER(iwp) ::  isign !< 
@@ -128,5 +129,4 @@
     INTEGER(iwp) ::  n     !< 
 
-    !  Array arguments 
     REAL(wp)     ::  a(*)     !< 
     REAL(wp)     ::  trigs(*) !< 
@@ -134,5 +134,4 @@
     INTEGER(iwp) ::  ifax(*)  !< 
 
-    !  Local scalars: 
     INTEGER(iwp) ::  i      !< 
     INTEGER(iwp) ::  ia     !< 
@@ -156,177 +155,185 @@
     INTEGER(iwp) ::  nx     !< 
 
-    !  Intrinsic functions 
-!    INTRINSIC MOD
-
-
-    !  Executable statements 
-
-    IF (ifax(10)/=n) CALL set99(trigs,ifax,n)
+
+    IF ( ifax(10) /= n )  CALL set99( trigs, ifax, n )
     nfax = ifax(1)
-    nx = n + 1
-    IF (MOD(n,2)==1) nx = n
-    nblox = 1 + (lot-1)/nfft
-    nvex = lot - (nblox-1)*nfft
-    IF (isign==-1) GO TO 50
-
-    ! isign=+1, spectral to gridpoint transform
-
-    istart = 1
-    DO nb = 1, nblox
-       ia = istart
-       i = istart
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO j = 1, nvex
-          a(i+inc) = 0.5_wp*a(i)
-          i = i + jump
-       END DO
-       IF (MOD(n,2)==1) GO TO 10
-       i = istart + n*inc
-       DO j = 1, nvex
-          a(i) = 0.5_wp*a(i)
-          i = i + jump
-       END DO
-10     CONTINUE
-       ia = istart + inc
-       la = 1
-       igo = + 1
-
-       DO k = 1, nfax
-          ifac = ifax(k+1)
-          ierr = -1
-          IF (igo==-1) GO TO 20
-          CALL rpassm(a(ia),a(ia+la*inc),work(1),work(ifac*la+1),trigs,inc,1, &
-               &          jump,nx,nvex,n,ifac,la,ierr)
-          GO TO 30
-20        CONTINUE
-          CALL rpassm(work(1),work(la+1),a(ia),a(ia+ifac*la*inc),trigs,1,inc,nx, &
-               &          jump,nvex,n,ifac,la,ierr)
-30        CONTINUE
-          IF (ierr/=0) GO TO 100
-          la = ifac*la
-          igo = -igo
+    nx   = n + 1
+    IF ( MOD(n,2) == 1 )  nx = n
+    nblox = 1 + (lot-1) / nfft
+    nvex = lot - (nblox-1) * nfft
+
+
+    IF ( isign == 1 )  THEN
+!
+!--    Backward fft: spectral to gridpoint transform
+
+       istart = 1
+       DO  nb = 1, nblox
+
           ia = istart
-       END DO
-
-       ! If necessary, copy results back to a
-
-       IF (MOD(nfax,2)==0) GO TO 40
-       ibase = 1
-       jbase = ia
-       DO jj = 1, nvex
-          i = ibase
-          j = jbase
-          DO ii = 1, n
-             a(j) = work(i)
-             i = i + 1
-             j = j + inc
-          END DO
-          ibase = ibase + nx
-          jbase = jbase + jump
-       END DO
-40     CONTINUE
-
-       ! Fill in zeros at end
-
-       ix = istart + n*inc
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO j = 1, nvex
-          a(ix) = 0.0_wp
-          a(ix+inc) = 0.0_wp
-          ix = ix + jump
-       END DO
-
-       istart = istart + nvex*jump
-       nvex = nfft
-    END DO
-    RETURN
-
-    ! isign=-1, gridpoint to spectral transform
-
-50  CONTINUE
-    istart = 1
-    DO nb = 1, nblox
-       ia = istart
-       la = n
-       igo = + 1
-
-       DO k = 1, nfax
-          ifac = ifax(nfax+2-k)
-          la = la/ifac
-          ierr = -1
-          IF (igo==-1) GO TO 60
-          CALL qpassm(a(ia),a(ia+ifac*la*inc),work(1),work(la+1),trigs,inc,1, &
-               &          jump,nx,nvex,n,ifac,la,ierr)
-          GO TO 70
-60        CONTINUE
-          CALL qpassm(work(1),work(ifac*la+1),a(ia),a(ia+la*inc),trigs,1,inc,nx, &
-               &          jump,nvex,n,ifac,la,ierr)
-70        CONTINUE
-          IF (ierr/=0) GO TO 100
-          igo = -igo
+          i = istart
+          !DIR$ IVDEP
+          !CDIR NODEP
+          !OCL NOVREC
+          DO  j = 1, nvex
+             a(i+inc) = 0.5_wp * a(i)
+             i = i + jump
+          ENDDO
+          IF ( MOD(n,2) /= 1 )  THEN
+             i = istart + n * inc
+             DO  j = 1, nvex
+                a(i) = 0.5_wp * a(i)
+                i = i + jump
+             ENDDO
+          ENDIF
           ia = istart + inc
-       END DO
-
-       ! If necessary, copy results back to a
-
-       IF (MOD(nfax,2)==0) GO TO 80
-       ibase = 1
-       jbase = ia
-       DO jj = 1, nvex
-          i = ibase
-          j = jbase
-          DO ii = 1, n
-             a(j) = work(i)
-             i = i + 1
-             j = j + inc
-          END DO
-          ibase = ibase + nx
-          jbase = jbase + jump
-       END DO
-80     CONTINUE
-
-       ! Shift a(0) & fill in zero imag parts
-
-       ix = istart
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO j = 1, nvex
-          a(ix) = a(ix+inc)
-          a(ix+inc) = 0.0_wp
-          ix = ix + jump
-       END DO
-       IF (MOD(n,2)==1) GO TO 90
-       iz = istart + (n+1)*inc
-       DO j = 1, nvex
-          a(iz) = 0.0_wp
-          iz = iz + jump
-       END DO
-90     CONTINUE
-
-       istart = istart + nvex*jump
-       nvex = nfft
-    END DO
-    RETURN
-
-    ! Error messages
-
-100 CONTINUE
-
-    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
-  END SUBROUTINE fft991cy
+          la = 1
+          igo = + 1
+
+          DO  k = 1, nfax
+
+             ifac = ifax(k+1)
+             ierr = -1
+             IF ( igo /= -1 )  THEN
+                CALL rpassm( a(ia), a(ia+la*inc), work(1), work(ifac*la+1), trigs, inc, 1, jump,   &
+                             nx, nvex, n, ifac, la, ierr )
+             ELSE
+                CALL rpassm( work(1), work(la+1), a(ia), a(ia+ifac*la*inc), trigs, 1, inc, nx,     &
+                             jump, nvex, 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
+             DO  jj = 1, nvex
+                i = ibase
+                j = jbase
+                DO  ii = 1, n
+                   a(j) = work(i)
+                   i = i + 1
+                   j = j + inc
+                ENDDO
+                ibase = ibase + nx
+                jbase = jbase + jump
+             ENDDO
+          ENDIF
+!
+!--       Fill in zeros at end
+          ix = istart + n*inc
+          !DIR$ IVDEP
+          !CDIR NODEP
+          !OCL NOVREC
+          DO  j = 1, nvex
+             a(ix) = 0.0_wp
+             a(ix+inc) = 0.0_wp
+             ix = ix + jump
+          ENDDO
+          istart = istart + nvex*jump
+          nvex = nfft
+
+       ENDDO
+
+    ELSEIF ( isign == -1 )  THEN
+!
+!--    Forward fft: gridpoint to spectral transform
+       istart = 1
+       DO  nb = 1, nblox
+          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( a(ia), a(ia+ifac*la*inc), work(1), work(la+1), trigs, inc, 1, jump,   &
+                             nx, nvex, n, ifac, la, ierr )
+             ELSE
+                CALL qpassm( work(1), work(ifac*la+1), a(ia), a(ia+la*inc), trigs, 1, inc, nx,     &
+                             jump, nvex, 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
+
+             igo = -igo
+             ia = istart + inc
+
+          ENDDO
+!
+!--       If necessary, copy results back to a
+          IF ( MOD(nfax,2) /= 0 )  THEN
+             ibase = 1
+             jbase = ia
+             DO  jj = 1, nvex
+                i = ibase
+                j = jbase
+                DO  ii = 1, n
+                   a(j) = work(i)
+                   i = i + 1
+                   j = j + inc
+                ENDDO
+                ibase = ibase + nx
+                jbase = jbase + jump
+             ENDDO
+          ENDIF
+!
+!--       Shift a(0) and fill in zero imag parts
+          ix = istart
+          !DIR$ IVDEP
+          !CDIR NODEP
+          !OCL NOVREC
+          DO  j = 1, nvex
+             a(ix) = a(ix+inc)
+             a(ix+inc) = 0.0_wp
+             ix = ix + jump
+          ENDDO
+          IF ( MOD(n,2) /= 1 )  THEN
+             iz = istart + (n+1) * inc
+             DO  j = 1, nvex
+                a(iz) = 0.0_wp
+                iz = iz + jump
+             ENDDO
+          ENDIF
+          istart = istart + nvex * jump
+          nvex = nfft
+
+       ENDDO
+
+    ENDIF
+
+ END SUBROUTINE fft991cy
 
 !------------------------------------------------------------------------------!
@@ -357,5 +364,5 @@
 !>         3 - ifac only catered for if la=n/ifac
 !------------------------------------------------------------------------------!
-  SUBROUTINE qpassm(a,b,c,d,trigs,inc1,inc2,inc3,inc4,lot,n,ifac,la,ierr)
+ SUBROUTINE qpassm( a, b, c, d, trigs, inc1, inc2, inc3, inc4, lot, n, ifac, la, ierr )
 
     USE kinds
@@ -363,5 +370,4 @@
     IMPLICIT NONE 
 
-    !  Scalar arguments 
     INTEGER(iwp) ::  ierr !< 
     INTEGER(iwp) ::  ifac !< 
@@ -374,6 +380,6 @@
     INTEGER(iwp) ::  n    !< 
 
-    !  Array arguments 
-    ! REAL :: a(n),b(n),c(n),d(n),trigs(n)
+!
+!-- Arrays are dimensioned with n
     REAL(wp) ::  a(*)     !< 
     REAL(wp) ::  b(*)     !< 
@@ -382,5 +388,4 @@
     REAL(wp) ::  trigs(*) !<
  
-    !  Local scalars: 
     REAL(wp) ::  a0     !< 
     REAL(wp) ::  a1     !< 
@@ -430,5 +435,4 @@
     INTEGER(iwp) ::  ia    !< 
     INTEGER(iwp) ::  ib    !< 
-    INTEGER(iwp) ::  ibad  !< 
     INTEGER(iwp) ::  ibase !< 
     INTEGER(iwp) ::  ic    !< 
@@ -461,779 +465,819 @@
     INTEGER(iwp) ::  m     !< 
 
-    !  Intrinsic functions 
-!    INTRINSIC REAL, SQRT
-
-    !  Data statements 
-    DATA sin36/0.587785252292473_wp/, sin72/0.951056516295154_wp/, &
-         &      qrt5/0.559016994374947_wp/, sin60/0.866025403784437_wp/
-
-
-    !  Executable statements 
-
-    m = n/ifac
-    iink = la*inc1
-    jink = la*inc2
-    ijump = (ifac-1)*iink
-    kstop = (n-ifac)/(2*ifac)
-
-    ibad = 1
-    IF (lot>nfft) GO TO 180
+
+    DATA  sin36/0.587785252292473_wp/, sin72/0.951056516295154_wp/,                                &
+          qrt5/0.559016994374947_wp/,  sin60/0.866025403784437_wp/
+
+
+    ierr = 0
+
+    IF ( lot > nfft )  THEN
+       ierr = 1
+       RETURN
+    ENDIF
+
+    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
-    ibad = 2
-    IF (igo<1 .OR. igo>6) GO TO 180
-    GO TO (10,40,70,100,130,160) igo
-
-    ! Coding for factor 2
-
-10  CONTINUE
-    ia = 1
-    ib = ia + iink
-    ja = 1
-    jb = ja + (2*m-la)*inc2
-
-    IF (la==m) GO TO 30
-
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = a(ia+i) + a(ib+i)
-          c(jb+j) = a(ia+i) - a(ib+i)
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-    ja = ja + jink
-    jink = 2*jink
-    jb = jb - jink
-    ibase = ibase + ijump
-    ijump = 2*ijump + iink
-    IF (ja==jb) GO TO 20
-    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
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-          DO ijk = 1, lot
-             c(ja+j) = a(ia+i) + (c1*a(ib+i)+s1*b(ib+i))
-             c(jb+j) = a(ia+i) - (c1*a(ib+i)+s1*b(ib+i))
-             d(ja+j) = (c1*b(ib+i)-s1*a(ib+i)) + b(ia+i)
-             d(jb+j) = (c1*b(ib+i)-s1*a(ib+i)) - b(ia+i)
-             i = i + inc3
-             j = j + inc4
-          END DO
-          ibase = ibase + inc1
-          jbase = jbase + inc2
-       END DO
-       ibase = ibase + ijump
-       ja = ja + jink
-       jb = jb - jink
-    END DO
-    IF (ja>jb) GO TO 170
-20  CONTINUE
-    jbase = 0
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = a(ia+i)
-          d(ja+j) = -a(ib+i)
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-
-    GO TO 170
-30  CONTINUE
-    z = 1.0_wp/REAL(n,KIND=wp)
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = z*(a(ia+i)+a(ib+i))
-          c(jb+j) = z*(a(ia+i)-a(ib+i))
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-    GO TO 170
-
-    ! Coding for factor 3
-
-40  CONTINUE
-    ia = 1
-    ib = ia + iink
-    ic = ib + iink
-    ja = 1
-    jb = ja + (2*m-la)*inc2
-    jc = jb
-
-    IF (la==m) GO TO 60
-
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
-          c(jb+j) = a(ia+i) - 0.5_wp*(a(ib+i)+a(ic+i))
-          d(jb+j) = sin60*(a(ic+i)-a(ib+i))
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-    ja = ja + jink
-    jink = 2*jink
-    jb = jb + jink
-    jc = jc - jink
-    ibase = ibase + ijump
-    ijump = 2*ijump + iink
-    IF (ja==jc) GO TO 50
-    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
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-          DO ijk = 1, lot
-             a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c2*a(ic+i)+s2*b(ic+i))
-             b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c2*b(ic+i)-s2*a(ic+i))
-             a2 = a(ia+i) - 0.5_wp*a1
-             b2 = b(ia+i) - 0.5_wp*b1
-             a3 = sin60*((c1*a(ib+i)+s1*b(ib+i))-(c2*a(ic+i)+s2*b(ic+i)))
-             b3 = sin60*((c1*b(ib+i)-s1*a(ib+i))-(c2*b(ic+i)-s2*a(ic+i)))
-             c(ja+j) = a(ia+i) + a1
-             d(ja+j) = b(ia+i) + b1
-             c(jb+j) = a2 + b3
-             d(jb+j) = b2 - a3
-             c(jc+j) = a2 - b3
-             d(jc+j) = -(b2+a3)
-             i = i + inc3
-             j = j + inc4
-          END DO
-          ibase = ibase + inc1
-          jbase = jbase + inc2
-       END DO
-       ibase = ibase + ijump
-       ja = ja + jink
-       jb = jb + jink
-       jc = jc - jink
-    END DO
-    IF (ja>jc) GO TO 170
-50  CONTINUE
-    jbase = 0
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = a(ia+i) + 0.5_wp*(a(ib+i)-a(ic+i))
-          d(ja+j) = -sin60*(a(ib+i)+a(ic+i))
-          c(jb+j) = a(ia+i) - (a(ib+i)-a(ic+i))
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-
-    GO TO 170
-60  CONTINUE
-    z = 1.0_wp/REAL(n,KIND=wp)
-    zsin60 = z*sin60
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = z*(a(ia+i)+(a(ib+i)+a(ic+i)))
-          c(jb+j) = z*(a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))
-          d(jb+j) = zsin60*(a(ic+i)-a(ib+i))
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-    GO TO 170
-
-    ! Coding for factor 4
-
-70  CONTINUE
-    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) GO TO 90
-
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))
-          c(jc+j) = (a(ia+i)+a(ic+i)) - (a(ib+i)+a(id+i))
-          c(jb+j) = a(ia+i) - a(ic+i)
-          d(jb+j) = a(id+i) - a(ib+i)
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-    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) GO TO 80
-    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
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-          DO ijk = 1, lot
-             a0 = a(ia+i) + (c2*a(ic+i)+s2*b(ic+i))
-             a2 = a(ia+i) - (c2*a(ic+i)+s2*b(ic+i))
-             a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c3*a(id+i)+s3*b(id+i))
-             a3 = (c1*a(ib+i)+s1*b(ib+i)) - (c3*a(id+i)+s3*b(id+i))
-             b0 = b(ia+i) + (c2*b(ic+i)-s2*a(ic+i))
-             b2 = b(ia+i) - (c2*b(ic+i)-s2*a(ic+i))
-             b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c3*b(id+i)-s3*a(id+i))
-             b3 = (c1*b(ib+i)-s1*a(ib+i)) - (c3*b(id+i)-s3*a(id+i))
-             c(ja+j) = a0 + a1
-             c(jc+j) = a0 - a1
-             d(ja+j) = b0 + b1
-             d(jc+j) = b1 - b0
-             c(jb+j) = a2 + b3
-             c(jd+j) = a2 - b3
-             d(jb+j) = b2 - a3
-             d(jd+j) = -(b2+a3)
-             i = i + inc3
-             j = j + inc4
-          END DO
-          ibase = ibase + inc1
-          jbase = jbase + inc2
-       END DO
-       ibase = ibase + ijump
-       ja = ja + jink
-       jb = jb + jink
-       jc = jc - jink
-       jd = jd - jink
-    END DO
-    IF (jb>jc) GO TO 170
-80  CONTINUE
-    sin45 = SQRT(0.5_wp)
-    jbase = 0
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = a(ia+i) + sin45*(a(ib+i)-a(id+i))
-          c(jb+j) = a(ia+i) - sin45*(a(ib+i)-a(id+i))
-          d(ja+j) = -a(ic+i) - sin45*(a(ib+i)+a(id+i))
-          d(jb+j) = a(ic+i) - sin45*(a(ib+i)+a(id+i))
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-
-    GO TO 170
-90  CONTINUE
-    z = 1.0_wp/REAL(n,KIND=wp)
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = z*((a(ia+i)+a(ic+i))+(a(ib+i)+a(id+i)))
-          c(jc+j) = z*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i)))
-          c(jb+j) = z*(a(ia+i)-a(ic+i))
-          d(jb+j) = z*(a(id+i)-a(ib+i))
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-    GO TO 170
-
-    ! Coding for factor 5
-
-100 CONTINUE
-    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) GO TO 120
-
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          a1 = a(ib+i) + a(ie+i)
-          a3 = a(ib+i) - a(ie+i)
-          a2 = a(ic+i) + a(id+i)
-          a4 = a(ic+i) - a(id+i)
-          a5 = a(ia+i) - 0.25_wp*(a1+a2)
-          a6 = qrt5*(a1-a2)
-          c(ja+j) = a(ia+i) + (a1+a2)
-          c(jb+j) = a5 + a6
-          c(jc+j) = a5 - a6
-          d(jb+j) = -sin72*a3 - sin36*a4
-          d(jc+j) = -sin36*a3 + sin72*a4
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-    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) GO TO 110
-    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
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-          DO ijk = 1, lot
-             a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c4*a(ie+i)+s4*b(ie+i))
-             a3 = (c1*a(ib+i)+s1*b(ib+i)) - (c4*a(ie+i)+s4*b(ie+i))
-             a2 = (c2*a(ic+i)+s2*b(ic+i)) + (c3*a(id+i)+s3*b(id+i))
-             a4 = (c2*a(ic+i)+s2*b(ic+i)) - (c3*a(id+i)+s3*b(id+i))
-             b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c4*b(ie+i)-s4*a(ie+i))
-             b3 = (c1*b(ib+i)-s1*a(ib+i)) - (c4*b(ie+i)-s4*a(ie+i))
-             b2 = (c2*b(ic+i)-s2*a(ic+i)) + (c3*b(id+i)-s3*a(id+i))
-             b4 = (c2*b(ic+i)-s2*a(ic+i)) - (c3*b(id+i)-s3*a(id+i))
-             a5 = a(ia+i) - 0.25_wp*(a1+a2)
-             a6 = qrt5*(a1-a2)
-             b5 = b(ia+i) - 0.25_wp*(b1+b2)
-             b6 = qrt5*(b1-b2)
-             a10 = a5 + a6
-             a20 = a5 - a6
-             b10 = b5 + b6
-             b20 = b5 - b6
-             a11 = sin72*b3 + sin36*b4
-             a21 = sin36*b3 - sin72*b4
-             b11 = sin72*a3 + sin36*a4
-             b21 = sin36*a3 - sin72*a4
-             c(ja+j) = a(ia+i) + (a1+a2)
-             c(jb+j) = a10 + a11
-             c(je+j) = a10 - a11
-             c(jc+j) = a20 + a21
-             c(jd+j) = a20 - a21
-             d(ja+j) = b(ia+i) + (b1+b2)
-             d(jb+j) = b10 - b11
-             d(je+j) = -(b10+b11)
-             d(jc+j) = b20 - b21
-             d(jd+j) = -(b20+b21)
-             i = i + inc3
-             j = j + inc4
-          END DO
-          ibase = ibase + inc1
-          jbase = jbase + inc2
-       END DO
-       ibase = ibase + ijump
-       ja = ja + jink
-       jb = jb + jink
-       jc = jc + jink
-       jd = jd - jink
-       je = je - jink
-    END DO
-    IF (jb>jd) GO TO 170
-110 CONTINUE
-    jbase = 0
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          a1 = a(ib+i) + a(ie+i)
-          a3 = a(ib+i) - a(ie+i)
-          a2 = a(ic+i) + a(id+i)
-          a4 = a(ic+i) - a(id+i)
-          a5 = a(ia+i) + 0.25_wp*(a3-a4)
-          a6 = qrt5*(a3+a4)
-          c(ja+j) = a5 + a6
-          c(jb+j) = a5 - a6
-          c(jc+j) = a(ia+i) - (a3-a4)
-          d(ja+j) = -sin36*a1 - sin72*a2
-          d(jb+j) = -sin72*a1 + sin36*a2
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-
-    GO TO 170
-120 CONTINUE
-    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
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          a1 = a(ib+i) + a(ie+i)
-          a3 = a(ib+i) - a(ie+i)
-          a2 = a(ic+i) + a(id+i)
-          a4 = a(ic+i) - a(id+i)
-          a5 = z*(a(ia+i)-0.25_wp*(a1+a2))
-          a6 = zqrt5*(a1-a2)
-          c(ja+j) = z*(a(ia+i)+(a1+a2))
-          c(jb+j) = a5 + a6
-          c(jc+j) = a5 - a6
-          d(jb+j) = -zsin72*a3 - zsin36*a4
-          d(jc+j) = -zsin36*a3 + zsin72*a4
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-    GO TO 170
-
-    ! Coding for factor 6
-
-130 CONTINUE
-    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) GO TO 150
-
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
-          c(ja+j) = (a(ia+i)+a(id+i)) + a11
-          c(jc+j) = (a(ia+i)+a(id+i)-0.5_wp*a11)
-          d(jc+j) = sin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
-          a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
-          c(jb+j) = (a(ia+i)-a(id+i)) - 0.5_wp*a11
-          d(jb+j) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
-          c(jd+j) = (a(ia+i)-a(id+i)) + a11
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-    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) GO TO 140
-    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
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-          DO ijk = 1, lot
-             a1 = c1*a(ib+i) + s1*b(ib+i)
-             b1 = c1*b(ib+i) - s1*a(ib+i)
-             a2 = c2*a(ic+i) + s2*b(ic+i)
-             b2 = c2*b(ic+i) - s2*a(ic+i)
-             a3 = c3*a(id+i) + s3*b(id+i)
-             b3 = c3*b(id+i) - s3*a(id+i)
-             a4 = c4*a(ie+i) + s4*b(ie+i)
-             b4 = c4*b(ie+i) - s4*a(ie+i)
-             a5 = c5*a(if+i) + s5*b(if+i)
-             b5 = c5*b(if+i) - s5*a(if+i)
-             a11 = (a2+a5) + (a1+a4)
-             a20 = (a(ia+i)+a3) - 0.5_wp*a11
-             a21 = sin60*((a2+a5)-(a1+a4))
-             b11 = (b2+b5) + (b1+b4)
-             b20 = (b(ia+i)+b3) - 0.5_wp*b11
-             b21 = sin60*((b2+b5)-(b1+b4))
-             c(ja+j) = (a(ia+i)+a3) + a11
-             d(ja+j) = (b(ia+i)+b3) + b11
-             c(jc+j) = a20 - b21
-             d(jc+j) = a21 + b20
-             c(je+j) = a20 + b21
-             d(je+j) = a21 - b20
-             a11 = (a2-a5) + (a4-a1)
-             a20 = (a(ia+i)-a3) - 0.5_wp*a11
-             a21 = sin60*((a4-a1)-(a2-a5))
-             b11 = (b5-b2) - (b4-b1)
-             b20 = (b3-b(ia+i)) - 0.5_wp*b11
-             b21 = sin60*((b5-b2)+(b4-b1))
-             c(jb+j) = a20 - b21
-             d(jb+j) = a21 - b20
-             c(jd+j) = a11 + (a(ia+i)-a3)
-             d(jd+j) = b11 + (b3-b(ia+i))
-             c(jf+j) = a20 + b21
-             d(jf+j) = a21 + b20
-             i = i + inc3
-             j = j + inc4
-          END DO
-          ibase = ibase + inc1
-          jbase = jbase + inc2
-       END DO
-       ibase = ibase + ijump
-       ja = ja + jink
-       jb = jb + jink
-       jc = jc + jink
-       jd = jd - jink
-       je = je - jink
-       jf = jf - jink
-    END DO
-    IF (jc>jd) GO TO 170
-140 CONTINUE
-    jbase = 0
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) + sin60*(a(ib+i)-a(if+i))
-          d(ja+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) - sin60*(a(ic+i)+a(ie+i))
-          c(jb+j) = a(ia+i) - (a(ic+i)-a(ie+i))
-          d(jb+j) = a(id+i) - (a(ib+i)+a(if+i))
-          c(jc+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) - sin60*(a(ib+i)-a(if+i))
-          d(jc+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) + sin60*(a(ic+i)+a(ie+i))
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-
-    GO TO 170
-150 CONTINUE
-    z = 1.0_wp/REAL(n,KIND=wp)
-    zsin60 = z*sin60
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
-          c(ja+j) = z*((a(ia+i)+a(id+i))+a11)
-          c(jc+j) = z*((a(ia+i)+a(id+i))-0.5_wp*a11)
-          d(jc+j) = zsin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
-          a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
-          c(jb+j) = z*((a(ia+i)-a(id+i))-0.5_wp*a11)
-          d(jb+j) = zsin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
-          c(jd+j) = z*((a(ia+i)-a(id+i))+a11)
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-    GO TO 170
-
-    ! Coding for factor 8
-
-160 CONTINUE
-    ibad = 3
-    IF (la/=m) GO TO 180
-    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
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = z*(((a(ia+i)+a(ie+i))+(a(ic+i)+a(ig+i)))+((a(id+i)+ &
-               &          a(ih+i))+(a(ib+i)+a(if+i))))
-          c(je+j) = z*(((a(ia+i)+a(ie+i))+(a(ic+i)+a(ig+i)))-((a(id+i)+ &
-               &          a(ih+i))+(a(ib+i)+a(if+i))))
-          c(jc+j) = z*((a(ia+i)+a(ie+i))-(a(ic+i)+a(ig+i)))
-          d(jc+j) = z*((a(id+i)+a(ih+i))-(a(ib+i)+a(if+i)))
-          c(jb+j) = z*(a(ia+i)-a(ie+i)) + zsin45*((a(ih+i)-a(id+i))-(a(if+ &
-               &          i)-a(ib+i)))
-          c(jd+j) = z*(a(ia+i)-a(ie+i)) - zsin45*((a(ih+i)-a(id+i))-(a(if+ &
-               &          i)-a(ib+i)))
-          d(jb+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) + &
-               &          z*(a(ig+i)-a(ic+i))
-          d(jd+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) - &
-               &          z*(a(ig+i)-a(ic+i))
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-
-    ! Return
-
-170 CONTINUE
-    ibad = 0
-180 CONTINUE
-    ierr = ibad
-    RETURN
-  END SUBROUTINE qpassm
+    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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = a(ia+i) + a(ib+i)
+                   c(jb+j) = a(ia+i) - a(ib+i)
+                   i = i + inc3
+                   j = j + inc4
+                ENDDO
+                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
+                      !DIR$ IVDEP
+                      !CDIR NODEP
+                      !OCL NOVREC
+                      DO  ijk = 1, lot
+                         c(ja+j) = a(ia+i) + (c1*a(ib+i)+s1*b(ib+i))
+                         c(jb+j) = a(ia+i) - (c1*a(ib+i)+s1*b(ib+i))
+                         d(ja+j) = (c1*b(ib+i)-s1*a(ib+i)) + b(ia+i)
+                         d(jb+j) = (c1*b(ib+i)-s1*a(ib+i)) - b(ia+i)
+                         i = i + inc3
+                         j = j + inc4
+                      ENDDO
+                      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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = a(ia+i)
+                   d(ja+j) = -a(ib+i)
+                   i = i + inc3
+                   j = j + inc4
+
+                ENDDO
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          ELSE
+
+             z = 1.0_wp/REAL(n,KIND=wp)
+             DO  l = 1, la
+                i = ibase
+                j = jbase
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = z*(a(ia+i)+a(ib+i))
+                   c(jb+j) = z*(a(ia+i)-a(ib+i))
+                   i = i + inc3
+                   j = j + inc4
+                ENDDO
+                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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
+                   c(jb+j) = a(ia+i) - 0.5_wp*(a(ib+i)+a(ic+i))
+                   d(jb+j) = sin60*(a(ic+i)-a(ib+i))
+                   i = i + inc3
+                   j = j + inc4
+                ENDDO
+                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
+                      !DIR$ IVDEP
+                      !CDIR NODEP
+                      !OCL NOVREC
+                      DO  ijk = 1, lot
+                         a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c2*a(ic+i)+s2*b(ic+i))
+                         b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c2*b(ic+i)-s2*a(ic+i))
+                         a2 = a(ia+i) - 0.5_wp*a1
+                         b2 = b(ia+i) - 0.5_wp*b1
+                         a3 = sin60*((c1*a(ib+i)+s1*b(ib+i))-(c2*a(ic+i)+s2*b(ic+i)))
+                         b3 = sin60*((c1*b(ib+i)-s1*a(ib+i))-(c2*b(ic+i)-s2*a(ic+i)))
+                         c(ja+j) = a(ia+i) + a1
+                         d(ja+j) = b(ia+i) + b1
+                         c(jb+j) = a2 + b3
+                         d(jb+j) = b2 - a3
+                         c(jc+j) = a2 - b3
+                         d(jc+j) = -(b2+a3)
+                         i = i + inc3
+                         j = j + inc4
+                      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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = a(ia+i) + 0.5_wp*(a(ib+i)-a(ic+i))
+                   d(ja+j) = -sin60*(a(ib+i)+a(ic+i))
+                   c(jb+j) = a(ia+i) - (a(ib+i)-a(ic+i))
+                   i = i + inc3
+                   j = j + inc4
+                ENDDO
+                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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = z*(a(ia+i)+(a(ib+i)+a(ic+i)))
+                   c(jb+j) = z*(a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))
+                   d(jb+j) = zsin60*(a(ic+i)-a(ib+i))
+                   i = i + inc3
+                   j = j + inc4
+                ENDDO
+                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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))
+                   c(jc+j) = (a(ia+i)+a(ic+i)) - (a(ib+i)+a(id+i))
+                   c(jb+j) = a(ia+i) - a(ic+i)
+                   d(jb+j) = a(id+i) - a(ib+i)
+                   i = i + inc3
+                   j = j + inc4
+                ENDDO
+                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
+                      !DIR$ IVDEP
+                      !CDIR NODEP
+                      !OCL NOVREC
+                      DO  ijk = 1, lot
+                         a0 = a(ia+i) + (c2*a(ic+i)+s2*b(ic+i))
+                         a2 = a(ia+i) - (c2*a(ic+i)+s2*b(ic+i))
+                         a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c3*a(id+i)+s3*b(id+i))
+                         a3 = (c1*a(ib+i)+s1*b(ib+i)) - (c3*a(id+i)+s3*b(id+i))
+                         b0 = b(ia+i) + (c2*b(ic+i)-s2*a(ic+i))
+                         b2 = b(ia+i) - (c2*b(ic+i)-s2*a(ic+i))
+                         b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c3*b(id+i)-s3*a(id+i))
+                         b3 = (c1*b(ib+i)-s1*a(ib+i)) - (c3*b(id+i)-s3*a(id+i))
+                         c(ja+j) = a0 + a1
+                         c(jc+j) = a0 - a1
+                         d(ja+j) = b0 + b1
+                         d(jc+j) = b1 - b0
+                         c(jb+j) = a2 + b3
+                         c(jd+j) = a2 - b3
+                         d(jb+j) = b2 - a3
+                         d(jd+j) = -(b2+a3)
+                         i = i + inc3
+                         j = j + inc4
+                      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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = a(ia+i) + sin45*(a(ib+i)-a(id+i))
+                   c(jb+j) = a(ia+i) - sin45*(a(ib+i)-a(id+i))
+                   d(ja+j) = -a(ic+i) - sin45*(a(ib+i)+a(id+i))
+                   d(jb+j) = a(ic+i) - sin45*(a(ib+i)+a(id+i))
+                   i = i + inc3
+                   j = j + inc4
+                ENDDO
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          ELSE
+
+             z = 1.0_wp / REAL( n, KIND=wp )
+             DO  l = 1, la
+                i = ibase
+                j = jbase
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = z*((a(ia+i)+a(ic+i))+(a(ib+i)+a(id+i)))
+                   c(jc+j) = z*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i)))
+                   c(jb+j) = z*(a(ia+i)-a(ic+i))
+                   d(jb+j) = z*(a(id+i)-a(ib+i))
+                   i = i + inc3
+                   j = j + inc4
+                ENDDO
+                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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   a1 = a(ib+i) + a(ie+i)
+                   a3 = a(ib+i) - a(ie+i)
+                   a2 = a(ic+i) + a(id+i)
+                   a4 = a(ic+i) - a(id+i)
+                   a5 = a(ia+i) - 0.25_wp*(a1+a2)
+                   a6 = qrt5*(a1-a2)
+                   c(ja+j) = a(ia+i) + (a1+a2)
+                   c(jb+j) = a5 + a6
+                   c(jc+j) = a5 - a6
+                   d(jb+j) = -sin72*a3 - sin36*a4
+                   d(jc+j) = -sin36*a3 + sin72*a4
+                   i = i + inc3
+                   j = j + inc4
+                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
+                      !DIR$ IVDEP
+                      !CDIR NODEP
+                      !OCL NOVREC
+                      DO  ijk = 1, lot
+                         a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c4*a(ie+i)+s4*b(ie+i))
+                         a3 = (c1*a(ib+i)+s1*b(ib+i)) - (c4*a(ie+i)+s4*b(ie+i))
+                         a2 = (c2*a(ic+i)+s2*b(ic+i)) + (c3*a(id+i)+s3*b(id+i))
+                         a4 = (c2*a(ic+i)+s2*b(ic+i)) - (c3*a(id+i)+s3*b(id+i))
+                         b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c4*b(ie+i)-s4*a(ie+i))
+                         b3 = (c1*b(ib+i)-s1*a(ib+i)) - (c4*b(ie+i)-s4*a(ie+i))
+                         b2 = (c2*b(ic+i)-s2*a(ic+i)) + (c3*b(id+i)-s3*a(id+i))
+                         b4 = (c2*b(ic+i)-s2*a(ic+i)) - (c3*b(id+i)-s3*a(id+i))
+                         a5 = a(ia+i) - 0.25_wp*(a1+a2)
+                         a6 = qrt5*(a1-a2)
+                         b5 = b(ia+i) - 0.25_wp*(b1+b2)
+                         b6 = qrt5*(b1-b2)
+                         a10 = a5 + a6
+                         a20 = a5 - a6
+                         b10 = b5 + b6
+                         b20 = b5 - b6
+                         a11 = sin72*b3 + sin36*b4
+                         a21 = sin36*b3 - sin72*b4
+                         b11 = sin72*a3 + sin36*a4
+                         b21 = sin36*a3 - sin72*a4
+                         c(ja+j) = a(ia+i) + (a1+a2)
+                         c(jb+j) = a10 + a11
+                         c(je+j) = a10 - a11
+                         c(jc+j) = a20 + a21
+                         c(jd+j) = a20 - a21
+                         d(ja+j) = b(ia+i) + (b1+b2)
+                         d(jb+j) = b10 - b11
+                         d(je+j) = -(b10+b11)
+                         d(jc+j) = b20 - b21
+                         d(jd+j) = -(b20+b21)
+                         i = i + inc3
+                         j = j + inc4
+                      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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   a1 = a(ib+i) + a(ie+i)
+                   a3 = a(ib+i) - a(ie+i)
+                   a2 = a(ic+i) + a(id+i)
+                   a4 = a(ic+i) - a(id+i)
+                   a5 = a(ia+i) + 0.25_wp*(a3-a4)
+                   a6 = qrt5*(a3+a4)
+                   c(ja+j) = a5 + a6
+                   c(jb+j) = a5 - a6
+                   c(jc+j) = a(ia+i) - (a3-a4)
+                   d(ja+j) = -sin36*a1 - sin72*a2
+                   d(jb+j) = -sin72*a1 + sin36*a2
+                   i = i + inc3
+                   j = j + inc4
+                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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   a1 = a(ib+i) + a(ie+i)
+                   a3 = a(ib+i) - a(ie+i)
+                   a2 = a(ic+i) + a(id+i)
+                   a4 = a(ic+i) - a(id+i)
+                   a5 = z*(a(ia+i)-0.25_wp*(a1+a2))
+                   a6 = zqrt5*(a1-a2)
+                   c(ja+j) = z*(a(ia+i)+(a1+a2))
+                   c(jb+j) = a5 + a6
+                   c(jc+j) = a5 - a6
+                   d(jb+j) = -zsin72*a3 - zsin36*a4
+                   d(jc+j) = -zsin36*a3 + zsin72*a4
+                   i = i + inc3
+                   j = j + inc4
+                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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
+                   c(ja+j) = (a(ia+i)+a(id+i)) + a11
+                   c(jc+j) = (a(ia+i)+a(id+i)-0.5_wp*a11)
+                   d(jc+j) = sin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
+                   a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
+                   c(jb+j) = (a(ia+i)-a(id+i)) - 0.5_wp*a11
+                   d(jb+j) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
+                   c(jd+j) = (a(ia+i)-a(id+i)) + a11
+                   i = i + inc3
+                   j = j + inc4
+                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
+             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
+                      !DIR$ IVDEP
+                      !CDIR NODEP
+                      !OCL NOVREC
+                      DO  ijk = 1, lot
+                         a1 = c1*a(ib+i) + s1*b(ib+i)
+                         b1 = c1*b(ib+i) - s1*a(ib+i)
+                         a2 = c2*a(ic+i) + s2*b(ic+i)
+                         b2 = c2*b(ic+i) - s2*a(ic+i)
+                         a3 = c3*a(id+i) + s3*b(id+i)
+                         b3 = c3*b(id+i) - s3*a(id+i)
+                         a4 = c4*a(ie+i) + s4*b(ie+i)
+                         b4 = c4*b(ie+i) - s4*a(ie+i)
+                         a5 = c5*a(if+i) + s5*b(if+i)
+                         b5 = c5*b(if+i) - s5*a(if+i)
+                         a11 = (a2+a5) + (a1+a4)
+                         a20 = (a(ia+i)+a3) - 0.5_wp*a11
+                         a21 = sin60*((a2+a5)-(a1+a4))
+                         b11 = (b2+b5) + (b1+b4)
+                         b20 = (b(ia+i)+b3) - 0.5_wp*b11
+                         b21 = sin60*((b2+b5)-(b1+b4))
+                         c(ja+j) = (a(ia+i)+a3) + a11
+                         d(ja+j) = (b(ia+i)+b3) + b11
+                         c(jc+j) = a20 - b21
+                         d(jc+j) = a21 + b20
+                         c(je+j) = a20 + b21
+                         d(je+j) = a21 - b20
+                         a11 = (a2-a5) + (a4-a1)
+                         a20 = (a(ia+i)-a3) - 0.5_wp*a11
+                         a21 = sin60*((a4-a1)-(a2-a5))
+                         b11 = (b5-b2) - (b4-b1)
+                         b20 = (b3-b(ia+i)) - 0.5_wp*b11
+                         b21 = sin60*((b5-b2)+(b4-b1))
+                         c(jb+j) = a20 - b21
+                         d(jb+j) = a21 - b20
+                         c(jd+j) = a11 + (a(ia+i)-a3)
+                         d(jd+j) = b11 + (b3-b(ia+i))
+                         c(jf+j) = a20 + b21
+                         d(jf+j) = a21 + b20
+                         i = i + inc3
+                         j = j + inc4
+                      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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) + sin60*(a(ib+i)-a(if+i))
+                   d(ja+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) - sin60*(a(ic+i)+a(ie+i))
+                   c(jb+j) = a(ia+i) - (a(ic+i)-a(ie+i))
+                   d(jb+j) = a(id+i) - (a(ib+i)+a(if+i))
+                   c(jc+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) - sin60*(a(ib+i)-a(if+i))
+                   d(jc+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) + sin60*(a(ic+i)+a(ie+i))
+                   i = i + inc3
+                   j = j + inc4
+                ENDDO
+                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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
+                   c(ja+j) = z*((a(ia+i)+a(id+i))+a11)
+                   c(jc+j) = z*((a(ia+i)+a(id+i))-0.5_wp*a11)
+                   d(jc+j) = zsin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
+                   a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
+                   c(jb+j) = z*((a(ia+i)-a(id+i))-0.5_wp*a11)
+                   d(jb+j) = zsin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
+                   c(jd+j) = z*((a(ia+i)-a(id+i))+a11)
+                   i = i + inc3
+                   j = j + inc4
+                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
+
+             !DIR$ IVDEP
+             !CDIR NODEP
+             !OCL NOVREC
+             DO  ijk = 1, lot
+                c(ja+j) = z*(((a(ia+i)+a(ie+i))+(a(ic+i)+a(ig+i)))+((a(id+i)+ a(ih+i))+(a(ib+i)+a(if+i))))
+                c(je+j) = z*(((a(ia+i)+a(ie+i))+(a(ic+i)+a(ig+i)))-((a(id+i)+ a(ih+i))+(a(ib+i)+a(if+i))))
+                c(jc+j) = z*((a(ia+i)+a(ie+i))-(a(ic+i)+a(ig+i)))
+                d(jc+j) = z*((a(id+i)+a(ih+i))-(a(ib+i)+a(if+i)))
+                c(jb+j) = z*(a(ia+i)-a(ie+i)) + zsin45*((a(ih+i)-a(id+i))-(a(if+i)-a(ib+i)))
+                c(jd+j) = z*(a(ia+i)-a(ie+i)) - zsin45*((a(ih+i)-a(id+i))-(a(if+i)-a(ib+i)))
+                d(jb+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) + z*(a(ig+i)-a(ic+i))
+                d(jd+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) - z*(a(ig+i)-a(ic+i))
+                i = i + inc3
+                j = j + inc4
+             ENDDO
+             ibase = ibase + inc1
+             jbase = jbase + inc2
+          ENDDO
+
+    END SELECT
+
+ END SUBROUTINE qpassm
 
 !------------------------------------------------------------------------------!
 ! Description:
 ! ------------
-!> @todo Missing subroutine description.
+!> Same as qpassm, but for backward fft
 !------------------------------------------------------------------------------!
-  SUBROUTINE rpassm(a,b,c,d,trigs,inc1,inc2,inc3,inc4,lot,n,ifac,la,ierr)
-    ! Dimension a(n),b(n),c(n),d(n),trigs(n)
+ SUBROUTINE rpassm( a, b, c, d, trigs, inc1, inc2, inc3, inc4, lot, n, ifac, la, ierr )
 
     USE kinds
@@ -1241,5 +1285,4 @@
     IMPLICIT NONE
 
-    !  Scalar arguments 
     INTEGER(iwp) ::  ierr !< 
     INTEGER(iwp) ::  ifac !< 
@@ -1251,6 +1294,6 @@
     INTEGER(iwp) ::  lot  !< 
     INTEGER(iwp) ::  n    !< 
-
-    !  Array arguments 
+!
+!-- Arrays are dimensioned with n
     REAL(wp) ::  a(*)     !< 
     REAL(wp) ::  b(*)     !< 
@@ -1259,5 +1302,4 @@
     REAL(wp) ::  trigs(*) !< 
 
-    !  Local scalars: 
     REAL(wp) ::  c1     !< 
     REAL(wp) ::  c2     !< 
@@ -1284,5 +1326,4 @@
     INTEGER(iwp) ::  ia    !< 
     INTEGER(iwp) ::  ib    !< 
-    INTEGER(iwp) ::  ibad  !< 
     INTEGER(iwp) ::  ibase !< 
     INTEGER(iwp) ::  ic    !< 
@@ -1315,5 +1356,4 @@
     INTEGER(iwp) ::  m     !< 
 
-    !  Local arrays: 
     REAL(wp) ::  a10(nfft) !< 
     REAL(wp) ::  a11(nfft) !< 
@@ -1325,758 +1365,770 @@
     REAL(wp) ::  b21(nfft) !< 
 
-    !  Intrinsic functions 
-!    INTRINSIC SQRT
-
-    !  Data statements 
-    DATA sin36/0.587785252292473_wp/, sin72/0.951056516295154_wp/, &
-         &      qrt5/0.559016994374947_wp/, sin60/0.866025403784437_wp/
-
-
-    !  Executable statements 
-
-    m = n/ifac
-    iink = la*inc1
-    jink = la*inc2
-    jump = (ifac-1)*jink
-    kstop = (n-ifac)/(2*ifac)
-
-    ibad = 1
-    IF (lot>nfft) GO TO 180
+
+    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)
+
+    IF ( lot > nfft )  THEN
+       ierr = 1
+       RETURN
+    ENDIF
     ibase = 0
     jbase = 0
     igo = ifac - 1
-    IF (igo==7) igo = 6
-    ibad = 2
-    IF (igo<1 .OR. igo>6) GO TO 180
-    GO TO (10,40,70,100,130,160) igo
-
-    ! Coding for factor 2
-
-10  CONTINUE
-    ia = 1
-    ib = ia + (2*m-la)*inc1
-    ja = 1
-    jb = ja + jink
-
-    IF (la==m) GO TO 30
-
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = a(ia+i) + a(ib+i)
-          c(jb+j) = a(ia+i) - a(ib+i)
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-    ia = ia + iink
-    iink = 2*iink
-    ib = ib - iink
-    ibase = 0
-    jbase = jbase + jump
-    jump = 2*jump + jink
-    IF (ia==ib) GO TO 20
-    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
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-          DO ijk = 1, lot
-             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))
-             i = i + inc3
-             j = j + inc4
-          END DO
-          ibase = ibase + inc1
-          jbase = jbase + inc2
-       END DO
-       ia = ia + iink
-       ib = ib - iink
-       jbase = jbase + jump
-    END DO
-    IF (ia>ib) GO TO 170
-20  CONTINUE
-    ibase = 0
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = a(ia+i)
-          c(jb+j) = -b(ia+i)
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-
-    GO TO 170
-30  CONTINUE
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i))
-          c(jb+j) = 2.0_wp*(a(ia+i)-a(ib+i))
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-    GO TO 170
-
-    ! Coding for factor 3
-
-40  CONTINUE
-    ia = 1
-    ib = ia + (2*m-la)*inc1
-    ic = ib
-    ja = 1
-    jb = ja + jink
-    jc = jb + jink
-
-    IF (la==m) GO TO 60
-
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = a(ia+i) + a(ib+i)
-          c(jb+j) = (a(ia+i)-0.5_wp*a(ib+i)) - (sin60*(b(ib+i)))
-          c(jc+j) = (a(ia+i)-0.5_wp*a(ib+i)) + (sin60*(b(ib+i)))
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-    ia = ia + iink
-    iink = 2*iink
-    ib = ib + iink
-    ic = ic - iink
-    jbase = jbase + jump
-    jump = 2*jump + jink
-    IF (ia==ic) GO TO 50
-    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
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-          DO ijk = 1, lot
-             c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
-             d(ja+j) = b(ia+i) + (b(ib+i)-b(ic+i))
-             c(jb+j) = c1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ &
-                  &            b(ic+i))))                                      &
-                  &    - s1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- &
-                  &            a(ic+i))))
-             d(jb+j) = s1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ &
-                  &            b(ic+i))))                                      &
-                  &    + c1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- &
-                  &            a(ic+i))))
-             c(jc+j) = c2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ &
-                  &            b(ic+i))))                                      &
-                  &    - s2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- &
-                  &            a(ic+i))))
-             d(jc+j) = s2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ &
-                  &            b(ic+i))))                                      &
-                  &    + c2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- &
-                  &            a(ic+i))))
-             i = i + inc3
-             j = j + inc4
-          END DO
-          ibase = ibase + inc1
-          jbase = jbase + inc2
-       END DO
-       ia = ia + iink
-       ib = ib + iink
-       ic = ic - iink
-       jbase = jbase + jump
-    END DO
-    IF (ia>ic) GO TO 170
-50  CONTINUE
-    ibase = 0
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = a(ia+i) + a(ib+i)
-          c(jb+j) = (0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
-          c(jc+j) = -(0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-
-    GO TO 170
-60  CONTINUE
-    ssin60 = 2.0_wp*sin60
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i))
-          c(jb+j) = (2.0_wp*a(ia+i)-a(ib+i)) - (ssin60*b(ib+i))
-          c(jc+j) = (2.0_wp*a(ia+i)-a(ib+i)) + (ssin60*b(ib+i))
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-    GO TO 170
-
-    ! Coding for factor 4
-
-70  CONTINUE
-    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) GO TO 90
-
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = (a(ia+i)+a(ic+i)) + a(ib+i)
-          c(jb+j) = (a(ia+i)-a(ic+i)) - b(ib+i)
-          c(jc+j) = (a(ia+i)+a(ic+i)) - a(ib+i)
-          c(jd+j) = (a(ia+i)-a(ic+i)) + b(ib+i)
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-    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) GO TO 80
-    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
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-          DO ijk = 1, lot
-             c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))
-             d(ja+j) = (b(ia+i)-b(ic+i)) + (b(ib+i)-b(id+i))
-             c(jc+j) = c2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) - s2*((b(ia+ &
-                  &            i)-b(ic+i))-(b(ib+i)-b(id+i)))
-             d(jc+j) = s2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) + c2*((b(ia+ &
-                  &            i)-b(ic+i))-(b(ib+i)-b(id+i)))
-             c(jb+j) = c1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) - s1*((b(ia+ &
-                  &            i)+b(ic+i))+(a(ib+i)-a(id+i)))
-             d(jb+j) = s1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) + c1*((b(ia+ &
-                  &            i)+b(ic+i))+(a(ib+i)-a(id+i)))
-             c(jd+j) = c3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) - s3*((b(ia+ &
-                  &            i)+b(ic+i))-(a(ib+i)-a(id+i)))
-             d(jd+j) = s3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) + c3*((b(ia+ &
-                  &            i)+b(ic+i))-(a(ib+i)-a(id+i)))
-             i = i + inc3
-             j = j + inc4
-          END DO
-          ibase = ibase + inc1
-          jbase = jbase + inc2
-       END DO
-       ia = ia + iink
-       ib = ib + iink
-       ic = ic - iink
-       id = id - iink
-       jbase = jbase + jump
-    END DO
-    IF (ib>ic) GO TO 170
-80  CONTINUE
-    ibase = 0
-    sin45 = SQRT(0.5_wp)
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = a(ia+i) + a(ib+i)
-          c(jb+j) = sin45*((a(ia+i)-a(ib+i))-(b(ia+i)+b(ib+i)))
-          c(jc+j) = b(ib+i) - b(ia+i)
-          c(jd+j) = -sin45*((a(ia+i)-a(ib+i))+(b(ia+i)+b(ib+i)))
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-
-    GO TO 170
-90  CONTINUE
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = 2.0_wp*((a(ia+i)+a(ic+i))+a(ib+i))
-          c(jb+j) = 2.0_wp*((a(ia+i)-a(ic+i))-b(ib+i))
-          c(jc+j) = 2.0_wp*((a(ia+i)+a(ic+i))-a(ib+i))
-          c(jd+j) = 2.0_wp*((a(ia+i)-a(ic+i))+b(ib+i))
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-
-    ! Coding for factor 5
-
-    GO TO 170
-100 CONTINUE
-    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) GO TO 120
-
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
-          c(jb+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) - &
-               &          (sin72*b(ib+i)+sin36*b(ic+i))
-          c(jc+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) - &
-               &          (sin36*b(ib+i)-sin72*b(ic+i))
-          c(jd+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) + &
-               &          (sin36*b(ib+i)-sin72*b(ic+i))
-          c(je+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) + &
-               &          (sin72*b(ib+i)+sin36*b(ic+i))
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-    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) GO TO 110
-    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
-!CDIR NODEP
-!OCL NOVREC
-          DO ijk = 1, lot
-
-             a10(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) + &
-                  &            qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
-             a20(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) - &
-                  &            qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
-             b10(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) + &
-                  &            qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
-             b20(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) - &
-                  &            qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
-             a11(ijk) = sin72*(b(ib+i)+b(ie+i)) + sin36*(b(ic+i)+b(id+i))
-             a21(ijk) = sin36*(b(ib+i)+b(ie+i)) - sin72*(b(ic+i)+b(id+i))
-             b11(ijk) = sin72*(a(ib+i)-a(ie+i)) + sin36*(a(ic+i)-a(id+i))
-             b21(ijk) = sin36*(a(ib+i)-a(ie+i)) - sin72*(a(ic+i)-a(id+i))
-
-             c(ja+j) = a(ia+i) + ((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))
-             d(ja+j) = b(ia+i) + ((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))
-             c(jb+j) = c1*(a10(ijk)-a11(ijk)) - s1*(b10(ijk)+b11(ijk))
-             d(jb+j) = s1*(a10(ijk)-a11(ijk)) + c1*(b10(ijk)+b11(ijk))
-             c(je+j) = c4*(a10(ijk)+a11(ijk)) - s4*(b10(ijk)-b11(ijk))
-             d(je+j) = s4*(a10(ijk)+a11(ijk)) + c4*(b10(ijk)-b11(ijk))
-             c(jc+j) = c2*(a20(ijk)-a21(ijk)) - s2*(b20(ijk)+b21(ijk))
-             d(jc+j) = s2*(a20(ijk)-a21(ijk)) + c2*(b20(ijk)+b21(ijk))
-             c(jd+j) = c3*(a20(ijk)+a21(ijk)) - s3*(b20(ijk)-b21(ijk))
-             d(jd+j) = s3*(a20(ijk)+a21(ijk)) + c3*(b20(ijk)-b21(ijk))
-
-             i = i + inc3
-             j = j + inc4
-          END DO
-          ibase = ibase + inc1
-          jbase = jbase + inc2
-       END DO
-       ia = ia + iink
-       ib = ib + iink
-       ic = ic + iink
-       id = id - iink
-       ie = ie - iink
-       jbase = jbase + jump
-    END DO
-    IF (ib>id) GO TO 170
-110 CONTINUE
-    ibase = 0
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = (a(ia+i)+a(ib+i)) + a(ic+i)
-          c(jb+j) = (qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - &
-               &          (sin36*b(ia+i)+sin72*b(ib+i))
-          c(je+j) = -(qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - &
-               &          (sin36*b(ia+i)+sin72*b(ib+i))
-          c(jc+j) = (qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - &
-               &          (sin72*b(ia+i)-sin36*b(ib+i))
-          c(jd+j) = -(qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - &
-               &          (sin72*b(ia+i)-sin36*b(ib+i))
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-
-    GO TO 170
-120 CONTINUE
-    qqrt5 = 2.0_wp*qrt5
-    ssin36 = 2.0_wp*sin36
-    ssin72 = 2.0_wp*sin72
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = 2.0_wp*(a(ia+i)+(a(ib+i)+a(ic+i)))
-          c(jb+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ &
-               &          i))) - (ssin72*b(ib+i)+ssin36*b(ic+i))
-          c(jc+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ &
-               &          i))) - (ssin36*b(ib+i)-ssin72*b(ic+i))
-          c(jd+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ &
-               &          i))) + (ssin36*b(ib+i)-ssin72*b(ic+i))
-          c(je+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ &
-               &          i))) + (ssin72*b(ib+i)+ssin36*b(ic+i))
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-    GO TO 170
-
-    ! Coding for factor 6
-
-130 CONTINUE
-    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) GO TO 150
-
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = (a(ia+i)+a(id+i)) + (a(ib+i)+a(ic+i))
-          c(jd+j) = (a(ia+i)-a(id+i)) - (a(ib+i)-a(ic+i))
-          c(jb+j) = ((a(ia+i)-a(id+i))+0.5_wp*(a(ib+i)-a(ic+i))) - (sin60*(b(ib+ &
-               &          i)+b(ic+i)))
-          c(jf+j) = ((a(ia+i)-a(id+i))+0.5_wp*(a(ib+i)-a(ic+i))) + (sin60*(b(ib+ &
-               &          i)+b(ic+i)))
-          c(jc+j) = ((a(ia+i)+a(id+i))-0.5_wp*(a(ib+i)+a(ic+i))) - (sin60*(b(ib+ &
-               &          i)-b(ic+i)))
-          c(je+j) = ((a(ia+i)+a(id+i))-0.5_wp*(a(ib+i)+a(ic+i))) + (sin60*(b(ib+ &
-               &          i)-b(ic+i)))
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-    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) GO TO 140
-    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
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-          DO ijk = 1, lot
-
-             a11(ijk) = (a(ie+i)+a(ib+i)) + (a(ic+i)+a(if+i))
-             a20(ijk) = (a(ia+i)+a(id+i)) - 0.5_wp*a11(ijk)
-             a21(ijk) = sin60*((a(ie+i)+a(ib+i))-(a(ic+i)+a(if+i)))
-             b11(ijk) = (b(ib+i)-b(ie+i)) + (b(ic+i)-b(if+i))
-             b20(ijk) = (b(ia+i)-b(id+i)) - 0.5_wp*b11(ijk)
-             b21(ijk) = sin60*((b(ib+i)-b(ie+i))-(b(ic+i)-b(if+i)))
-
-             c(ja+j) = (a(ia+i)+a(id+i)) + a11(ijk)
-             d(ja+j) = (b(ia+i)-b(id+i)) + b11(ijk)
-             c(jc+j) = c2*(a20(ijk)-b21(ijk)) - s2*(b20(ijk)+a21(ijk))
-             d(jc+j) = s2*(a20(ijk)-b21(ijk)) + c2*(b20(ijk)+a21(ijk))
-             c(je+j) = c4*(a20(ijk)+b21(ijk)) - s4*(b20(ijk)-a21(ijk))
-             d(je+j) = s4*(a20(ijk)+b21(ijk)) + c4*(b20(ijk)-a21(ijk))
-
-             a11(ijk) = (a(ie+i)-a(ib+i)) + (a(ic+i)-a(if+i))
-             b11(ijk) = (b(ie+i)+b(ib+i)) - (b(ic+i)+b(if+i))
-             a20(ijk) = (a(ia+i)-a(id+i)) - 0.5_wp*a11(ijk)
-             a21(ijk) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
-             b20(ijk) = (b(ia+i)+b(id+i)) + 0.5_wp*b11(ijk)
-             b21(ijk) = sin60*((b(ie+i)+b(ib+i))+(b(ic+i)+b(if+i)))
-
-             c(jd+j) = c3*((a(ia+i)-a(id+i))+a11(ijk)) - s3*((b(ia+i)+b(id+ &
-                  &            i))-b11(ijk))
-             d(jd+j) = s3*((a(ia+i)-a(id+i))+a11(ijk)) + c3*((b(ia+i)+b(id+ &
-                  &            i))-b11(ijk))
-             c(jb+j) = c1*(a20(ijk)-b21(ijk)) - s1*(b20(ijk)-a21(ijk))
-             d(jb+j) = s1*(a20(ijk)-b21(ijk)) + c1*(b20(ijk)-a21(ijk))
-             c(jf+j) = c5*(a20(ijk)+b21(ijk)) - s5*(b20(ijk)+a21(ijk))
-             d(jf+j) = s5*(a20(ijk)+b21(ijk)) + c5*(b20(ijk)+a21(ijk))
-
-             i = i + inc3
-             j = j + inc4
-          END DO
-          ibase = ibase + inc1
-          jbase = jbase + inc2
-       END DO
-       ia = ia + iink
-       ib = ib + iink
-       ic = ic + iink
-       id = id - iink
-       ie = ie - iink
-       if = if - iink
-       jbase = jbase + jump
-    END DO
-    IF (ic>id) GO TO 170
-140 CONTINUE
-    ibase = 0
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = a(ib+i) + (a(ia+i)+a(ic+i))
-          c(jd+j) = b(ib+i) - (b(ia+i)+b(ic+i))
-          c(jb+j) = (sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i))
-          c(jf+j) = -(sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i))
-          c(jc+j) = sin60*(b(ic+i)-b(ia+i)) + (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i))
-          c(je+j) = sin60*(b(ic+i)-b(ia+i)) - (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i))
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-
-    GO TO 170
-150 CONTINUE
-    ssin60 = 2.0_wp*sin60
-    DO l = 1, la
-       i = ibase
-       j = jbase
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = (2.0_wp*(a(ia+i)+a(id+i))) + (2.0_wp*(a(ib+i)+a(ic+i)))
-          c(jd+j) = (2.0_wp*(a(ia+i)-a(id+i))) - (2.0_wp*(a(ib+i)-a(ic+i)))
-          c(jb+j) = (2.0_wp*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) - (ssin60*(b(ib+ &
-               &          i)+b(ic+i)))
-          c(jf+j) = (2.0_wp*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) + (ssin60*(b(ib+ &
-               &          i)+b(ic+i)))
-          c(jc+j) = (2.0_wp*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) - (ssin60*(b(ib+ &
-               &          i)-b(ic+i)))
-          c(je+j) = (2.0_wp*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) + (ssin60*(b(ib+ &
-               &          i)-b(ic+i)))
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-    GO TO 170
-
-    ! Coding for factor 8
-
-160 CONTINUE
-    ibad = 3
-    IF (la/=m) GO TO 180
-    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
-!DIR$ IVDEP
-!CDIR NODEP
-!OCL NOVREC
-       DO ijk = 1, lot
-          c(ja+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))+(a(ib+i)+a(id+i)))
-          c(je+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))-(a(ib+i)+a(id+i)))
-          c(jc+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))-(b(ib+i)-b(id+i)))
-          c(jg+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))+(b(ib+i)-b(id+i)))
-          c(jb+j) = 2.0_wp*((a(ia+i)-a(ie+i))-b(ic+i)) + ssin45*((a(ib+i)-a(id+ &
-               &              i))-(b(ib+i)+b(id+i)))
-          c(jf+j) = 2.0_wp*((a(ia+i)-a(ie+i))-b(ic+i)) - ssin45*((a(ib+i)-a(id+ &
-               &              i))-(b(ib+i)+b(id+i)))
-          c(jd+j) = 2.0_wp*((a(ia+i)-a(ie+i))+b(ic+i)) - ssin45*((a(ib+i)-a(id+ &
-               &              i))+(b(ib+i)+b(id+i)))
-          c(jh+j) = 2.0_wp*((a(ia+i)-a(ie+i))+b(ic+i)) + ssin45*((a(ib+i)-a(id+ &
-               &              i))+(b(ib+i)+b(id+i)))
-          i = i + inc3
-          j = j + inc4
-       END DO
-       ibase = ibase + inc1
-       jbase = jbase + inc2
-    END DO
-
-    ! Return
-
-170 CONTINUE
-    ibad = 0
-180 CONTINUE
-    ierr = ibad
-    RETURN
-  END SUBROUTINE rpassm
+    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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = a(ia+i) + a(ib+i)
+                   c(jb+j) = a(ia+i) - a(ib+i)
+                   i = i + inc3
+                   j = j + inc4
+                ENDDO
+                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
+                      !DIR$ IVDEP
+                      !CDIR NODEP
+                      !OCL NOVREC
+                      DO  ijk = 1, lot
+                         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))
+                         i = i + inc3
+                         j = j + inc4
+                      ENDDO
+                      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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = a(ia+i)
+                   c(jb+j) = -b(ia+i)
+                   i = i + inc3
+                   j = j + inc4
+                ENDDO
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          ELSE
+
+             DO  l = 1, la
+                i = ibase
+                j = jbase
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i))
+                   c(jb+j) = 2.0_wp*(a(ia+i)-a(ib+i))
+                   i = i + inc3
+                   j = j + inc4
+                ENDDO
+                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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = a(ia+i) + a(ib+i)
+                   c(jb+j) = (a(ia+i)-0.5_wp*a(ib+i)) - (sin60*(b(ib+i)))
+                   c(jc+j) = (a(ia+i)-0.5_wp*a(ib+i)) + (sin60*(b(ib+i)))
+                   i = i + inc3
+                   j = j + inc4
+                ENDDO
+                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
+                      !DIR$ IVDEP
+                      !CDIR NODEP
+                      !OCL NOVREC
+                      DO  ijk = 1, lot
+                         c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
+                         d(ja+j) = b(ia+i) + (b(ib+i)-b(ic+i))
+                         c(jb+j) = c1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ b(ic+i)))) &
+                                   - s1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- a(ic+i))))
+                         d(jb+j) = s1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ b(ic+i)))) &
+                                   + c1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- a(ic+i))))
+                         c(jc+j) = c2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ b(ic+i)))) &
+                                   - s2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- a(ic+i))))
+                         d(jc+j) = s2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ b(ic+i)))) &
+                                   + c2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- a(ic+i))))
+                         i = i + inc3
+                         j = j + inc4
+                      ENDDO
+                      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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = a(ia+i) + a(ib+i)
+                   c(jb+j) = (0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
+                   c(jc+j) = -(0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
+                   i = i + inc3
+                   j = j + inc4
+                ENDDO
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          ELSE
+
+             ssin60 = 2.0_wp * sin60
+             DO  l = 1, la
+                i = ibase
+                j = jbase
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i))
+                   c(jb+j) = (2.0_wp*a(ia+i)-a(ib+i)) - (ssin60*b(ib+i))
+                   c(jc+j) = (2.0_wp*a(ia+i)-a(ib+i)) + (ssin60*b(ib+i))
+                   i = i + inc3
+                   j = j + inc4
+                ENDDO
+                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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = (a(ia+i)+a(ic+i)) + a(ib+i)
+                   c(jb+j) = (a(ia+i)-a(ic+i)) - b(ib+i)
+                   c(jc+j) = (a(ia+i)+a(ic+i)) - a(ib+i)
+                   c(jd+j) = (a(ia+i)-a(ic+i)) + b(ib+i)
+                   i = i + inc3
+                   j = j + inc4
+               ENDDO
+                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
+                      !DIR$ IVDEP
+                      !CDIR NODEP
+                      !OCL NOVREC
+                      DO  ijk = 1, lot
+                         c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))
+                         d(ja+j) = (b(ia+i)-b(ic+i)) + (b(ib+i)-b(id+i))
+                         c(jc+j) = c2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) - s2*((b(ia+i)-b(ic+i))&
+                                   -(b(ib+i)-b(id+i)))
+                         d(jc+j) = s2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) + c2*((b(ia+i)-b(ic+i))&
+                                   -(b(ib+i)-b(id+i)))
+                         c(jb+j) = c1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) - s1*((b(ia+i)+b(ic+i))&
+                                   +(a(ib+i)-a(id+i)))
+                         d(jb+j) = s1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) + c1*((b(ia+i)+b(ic+i))&
+                                   +(a(ib+i)-a(id+i)))
+                         c(jd+j) = c3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) - s3*((b(ia+i)+b(ic+i))&
+                                   -(a(ib+i)-a(id+i)))
+                         d(jd+j) = s3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) + c3*((b(ia+i)+b(ic+i))&
+                                   -(a(ib+i)-a(id+i)))
+                         i = i + inc3
+                         j = j + inc4
+                      ENDDO
+                      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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = a(ia+i) + a(ib+i)
+                   c(jb+j) = sin45*((a(ia+i)-a(ib+i))-(b(ia+i)+b(ib+i)))
+                   c(jc+j) = b(ib+i) - b(ia+i)
+                   c(jd+j) = -sin45*((a(ia+i)-a(ib+i))+(b(ia+i)+b(ib+i)))
+                   i = i + inc3
+                   j = j + inc4
+                ENDDO
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          ELSE
+
+             DO  l = 1, la
+                i = ibase
+                j = jbase
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = 2.0_wp*((a(ia+i)+a(ic+i))+a(ib+i))
+                   c(jb+j) = 2.0_wp*((a(ia+i)-a(ic+i))-b(ib+i))
+                   c(jc+j) = 2.0_wp*((a(ia+i)+a(ic+i))-a(ib+i))
+                   c(jd+j) = 2.0_wp*((a(ia+i)-a(ic+i))+b(ib+i))
+                   i = i + inc3
+                   j = j + inc4
+                ENDDO
+                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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
+                   c(jb+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) -        &
+                             (sin72*b(ib+i)+sin36*b(ic+i))
+                   c(jc+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) -        &
+                             (sin36*b(ib+i)-sin72*b(ic+i))
+                   c(jd+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) +        &
+                             (sin36*b(ib+i)-sin72*b(ic+i))
+                   c(je+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) +        &
+                             (sin72*b(ib+i)+sin36*b(ic+i))
+                   i = i + inc3
+                   j = j + inc4
+                ENDDO
+                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
+                      !CDIR NODEP
+                      !OCL NOVREC
+                      DO  ijk = 1, lot
+                         a10(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) +      &
+                                    qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
+                         a20(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) -      &
+                                    qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
+                         b10(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) +      &
+                                    qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
+                         b20(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) -      &
+                                    qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
+                         a11(ijk) = sin72*(b(ib+i)+b(ie+i)) + sin36*(b(ic+i)+b(id+i))
+                         a21(ijk) = sin36*(b(ib+i)+b(ie+i)) - sin72*(b(ic+i)+b(id+i))
+                         b11(ijk) = sin72*(a(ib+i)-a(ie+i)) + sin36*(a(ic+i)-a(id+i))
+                         b21(ijk) = sin36*(a(ib+i)-a(ie+i)) - sin72*(a(ic+i)-a(id+i))
+                         c(ja+j) = a(ia+i) + ((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))
+                         d(ja+j) = b(ia+i) + ((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))
+                         c(jb+j) = c1*(a10(ijk)-a11(ijk)) - s1*(b10(ijk)+b11(ijk))
+                         d(jb+j) = s1*(a10(ijk)-a11(ijk)) + c1*(b10(ijk)+b11(ijk))
+                         c(je+j) = c4*(a10(ijk)+a11(ijk)) - s4*(b10(ijk)-b11(ijk))
+                         d(je+j) = s4*(a10(ijk)+a11(ijk)) + c4*(b10(ijk)-b11(ijk))
+                         c(jc+j) = c2*(a20(ijk)-a21(ijk)) - s2*(b20(ijk)+b21(ijk))
+                         d(jc+j) = s2*(a20(ijk)-a21(ijk)) + c2*(b20(ijk)+b21(ijk))
+                         c(jd+j) = c3*(a20(ijk)+a21(ijk)) - s3*(b20(ijk)-b21(ijk))
+                         d(jd+j) = s3*(a20(ijk)+a21(ijk)) + c3*(b20(ijk)-b21(ijk))
+
+                         i = i + inc3
+                         j = j + inc4
+                      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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = (a(ia+i)+a(ib+i)) + a(ic+i)
+                   c(jb+j) = (qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) -        &
+                             (sin36*b(ia+i)+sin72*b(ib+i))
+                   c(je+j) = -(qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) -       &
+                             (sin36*b(ia+i)+sin72*b(ib+i))
+                   c(jc+j) = (qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) -        &
+                             (sin72*b(ia+i)-sin36*b(ib+i))
+                   c(jd+j) = -(qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) -       &
+                             (sin72*b(ia+i)-sin36*b(ib+i))
+                   i = i + inc3
+                   j = j + inc4
+                ENDDO
+                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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = 2.0_wp*(a(ia+i)+(a(ib+i)+a(ic+i)))
+                   c(jb+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+i)))  &
+                             - (ssin72*b(ib+i)+ssin36*b(ic+i))
+                   c(jc+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+i)))  &
+                             - (ssin36*b(ib+i)-ssin72*b(ic+i))
+                   c(jd+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+i)))  &
+                             + (ssin36*b(ib+i)-ssin72*b(ic+i))
+                   c(je+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+i)))  &
+                             + (ssin72*b(ib+i)+ssin36*b(ic+i))
+                   i = i + inc3
+                   j = j + inc4
+                ENDDO
+                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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = (a(ia+i)+a(id+i)) + (a(ib+i)+a(ic+i))
+                   c(jd+j) = (a(ia+i)-a(id+i)) - (a(ib+i)-a(ic+i))
+                   c(jb+j) = ((a(ia+i)-a(id+i))+0.5_wp*(a(ib+i)-a(ic+i))) - (sin60*(b(ib+i)+b(ic+i)))
+                   c(jf+j) = ((a(ia+i)-a(id+i))+0.5_wp*(a(ib+i)-a(ic+i))) + (sin60*(b(ib+i)+b(ic+i)))
+                   c(jc+j) = ((a(ia+i)+a(id+i))-0.5_wp*(a(ib+i)+a(ic+i))) - (sin60*(b(ib+i)-b(ic+i)))
+                   c(je+j) = ((a(ia+i)+a(id+i))-0.5_wp*(a(ib+i)+a(ic+i))) + (sin60*(b(ib+i)-b(ic+i)))
+                   i = i + inc3
+                   j = j + inc4
+                ENDDO
+                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
+                      !DIR$ IVDEP
+                      !CDIR NODEP
+                      !OCL NOVREC
+                      DO  ijk = 1, lot
+                         a11(ijk) = (a(ie+i)+a(ib+i)) + (a(ic+i)+a(if+i))
+                         a20(ijk) = (a(ia+i)+a(id+i)) - 0.5_wp*a11(ijk)
+                         a21(ijk) = sin60*((a(ie+i)+a(ib+i))-(a(ic+i)+a(if+i)))
+                         b11(ijk) = (b(ib+i)-b(ie+i)) + (b(ic+i)-b(if+i))
+                         b20(ijk) = (b(ia+i)-b(id+i)) - 0.5_wp*b11(ijk)
+                         b21(ijk) = sin60*((b(ib+i)-b(ie+i))-(b(ic+i)-b(if+i)))
+
+                         c(ja+j) = (a(ia+i)+a(id+i)) + a11(ijk)
+                         d(ja+j) = (b(ia+i)-b(id+i)) + b11(ijk)
+                         c(jc+j) = c2*(a20(ijk)-b21(ijk)) - s2*(b20(ijk)+a21(ijk))
+                         d(jc+j) = s2*(a20(ijk)-b21(ijk)) + c2*(b20(ijk)+a21(ijk))
+                         c(je+j) = c4*(a20(ijk)+b21(ijk)) - s4*(b20(ijk)-a21(ijk))
+                         d(je+j) = s4*(a20(ijk)+b21(ijk)) + c4*(b20(ijk)-a21(ijk))
+
+                         a11(ijk) = (a(ie+i)-a(ib+i)) + (a(ic+i)-a(if+i))
+                         b11(ijk) = (b(ie+i)+b(ib+i)) - (b(ic+i)+b(if+i))
+                         a20(ijk) = (a(ia+i)-a(id+i)) - 0.5_wp*a11(ijk)
+                         a21(ijk) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
+                         b20(ijk) = (b(ia+i)+b(id+i)) + 0.5_wp*b11(ijk)
+                         b21(ijk) = sin60*((b(ie+i)+b(ib+i))+(b(ic+i)+b(if+i)))
+
+                         c(jd+j) = c3*((a(ia+i)-a(id+i))+a11(ijk)) - s3*((b(ia+i)+b(id+i))-b11(ijk))
+                         d(jd+j) = s3*((a(ia+i)-a(id+i))+a11(ijk)) + c3*((b(ia+i)+b(id+i))-b11(ijk))
+                         c(jb+j) = c1*(a20(ijk)-b21(ijk)) - s1*(b20(ijk)-a21(ijk))
+                         d(jb+j) = s1*(a20(ijk)-b21(ijk)) + c1*(b20(ijk)-a21(ijk))
+                         c(jf+j) = c5*(a20(ijk)+b21(ijk)) - s5*(b20(ijk)+a21(ijk))
+                         d(jf+j) = s5*(a20(ijk)+b21(ijk)) + c5*(b20(ijk)+a21(ijk))
+
+                         i = i + inc3
+                         j = j + inc4
+                      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
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = a(ib+i) + (a(ia+i)+a(ic+i))
+                   c(jd+j) = b(ib+i) - (b(ia+i)+b(ic+i))
+                   c(jb+j) = (sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i))
+                   c(jf+j) = -(sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i))
+                   c(jc+j) = sin60*(b(ic+i)-b(ia+i)) + (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i))
+                   c(je+j) = sin60*(b(ic+i)-b(ia+i)) - (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i))
+                   i = i + inc3
+                   j = j + inc4
+                ENDDO
+                ibase = ibase + inc1
+                jbase = jbase + inc2
+             ENDDO
+
+          ELSE
+
+             ssin60 = 2.0_wp * sin60
+             DO  l = 1, la
+                i = ibase
+                j = jbase
+                !DIR$ IVDEP
+                !CDIR NODEP
+                !OCL NOVREC
+                DO  ijk = 1, lot
+                   c(ja+j) = (2.0_wp*(a(ia+i)+a(id+i))) + (2.0_wp*(a(ib+i)+a(ic+i)))
+                   c(jd+j) = (2.0_wp*(a(ia+i)-a(id+i))) - (2.0_wp*(a(ib+i)-a(ic+i)))
+                   c(jb+j) = (2.0_wp*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) - (ssin60*(b(ib+i)+b(ic+i)))
+                   c(jf+j) = (2.0_wp*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) + (ssin60*(b(ib+i)+b(ic+i)))
+                   c(jc+j) = (2.0_wp*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) - (ssin60*(b(ib+i)-b(ic+i)))
+                   c(je+j) = (2.0_wp*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) + (ssin60*(b(ib+i)-b(ic+i)))
+                   i = i + inc3
+                   j = j + inc4
+                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 + 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
+             !DIR$ IVDEP
+             !CDIR NODEP
+             !OCL NOVREC
+             DO  ijk = 1, lot
+                c(ja+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))+(a(ib+i)+a(id+i)))
+                c(je+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))-(a(ib+i)+a(id+i)))
+                c(jc+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))-(b(ib+i)-b(id+i)))
+                c(jg+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))+(b(ib+i)-b(id+i)))
+                c(jb+j) = 2.0_wp*((a(ia+i)-a(ie+i))-b(ic+i)) + ssin45*((a(ib+i)-a(id+i))-(b(ib+i)+b(id+i)))
+                c(jf+j) = 2.0_wp*((a(ia+i)-a(ie+i))-b(ic+i)) - ssin45*((a(ib+i)-a(id+i))-(b(ib+i)+b(id+i)))
+                c(jd+j) = 2.0_wp*((a(ia+i)-a(ie+i))+b(ic+i)) - ssin45*((a(ib+i)-a(id+i))+(b(ib+i)+b(id+i)))
+                c(jh+j) = 2.0_wp*((a(ia+i)-a(ie+i))+b(ic+i)) + ssin45*((a(ib+i)-a(id+i))+(b(ib+i)+b(id+i)))
+                i = i + inc3
+                j = j + inc4
+             ENDDO
+             ibase = ibase + inc1
+             jbase = jbase + inc2
+          ENDDO
+
+    END SELECT
+
+ END SUBROUTINE rpassm
 
 !------------------------------------------------------------------------------!
@@ -2088,8 +2140,8 @@
 !> functins required by fft99 & fft991cy
 !------------------------------------------------------------------------------!
-  SUBROUTINE set99(trigs,ifax,n)
-
-
-    USE control_parameters,                                                    &
+ SUBROUTINE set99( trigs, ifax, n )
+
+
+    USE control_parameters,                                                                        &
         ONLY:  message_string
 
@@ -2098,13 +2150,10 @@
     IMPLICIT NONE
 
-    !  Scalar arguments 
     INTEGER(iwp) ::  n !< 
 
-    !  Array arguments 
     INTEGER(iwp) ::  ifax(*)  !< 
     REAL(wp)     ::  trigs(*) !< 
 
 
-    !  Local scalars: 
     REAL(wp) ::  angle    !< 
     REAL(wp) ::  del      !< 
@@ -2119,68 +2168,64 @@
     INTEGER(iwp) ::  nu   !< 
 
-    !  Local arrays: 
     INTEGER(iwp) ::  jfax(10) !< 
     INTEGER(iwp) ::  lfax(7)  !< 
 
-    !  Intrinsic functions 
-!    INTRINSIC ASIN, COS, MOD, REAL, SIN
-
-    !  Data statements 
-    DATA lfax/6, 8, 5, 4, 3, 2, 1/
-
-
-    !  Executable statements 
+    DATA  lfax / 6, 8, 5, 4, 3, 2, 1 /
+
+
     ixxx = 1
 
-    del = 4.0_wp*ASIN(1.0_wp)/REAL(n,KIND=wp)
+    del = 4.0_wp * ASIN( 1.0_wp ) / REAL( n, KIND=wp )
     nil = 0
     nhl = (n/2) - 1
-    DO k = nil, nhl
-       angle = REAL(k,KIND=wp)*del
-       trigs(2*k+1) = COS(angle)
-       trigs(2*k+2) = SIN(angle)
-    END DO
-
-    ! Find factors of n (8,6,5,4,3,2; only one 8 allowed)
-    ! Look for sixes first, store factors in descending order
-    nu = n
+    DO  k = nil, nhl
+       angle = REAL( k, KIND=wp ) * del
+       trigs(2*k+1) = COS( angle )
+       trigs(2*k+2) = SIN( angle )
+    ENDDO
+
+!
+!-- Find factors of n (8,6,5,4,3,2; only one 8 allowed).
+!-- Look for sixes first, store factors in descending order.
+    nu   = n
     ifac = 6
-    k = 0
-    l = 1
-10  CONTINUE
-    IF (MOD(nu,ifac)/=0) GO TO 30
-    k = k + 1
-    jfax(k) = ifac
-    IF (ifac/=8) GO TO 20
-    IF (k==1) GO TO 20
-    jfax(1) = 8
-    jfax(k) = 6
-20  CONTINUE
-    nu = nu/ifac
-    IF (nu==1) GO TO 40
-    IF (ifac/=8) GO TO 10
-30  CONTINUE
-    l = l + 1
-    ifac = lfax(l)
-    IF (ifac>1) GO TO 10
-
-!    WRITE (nout,'(A,I4,A)') ' n =',n,' - Contains illegal factors'
-    message_string = 'number of gridpoints along x or/and y ' //               &
-                     'contain illegal  factors' //                             &
-                     '&only factors 2, 3, 5 are allowed' 
-    CALL message( 'temperton_fft', 'PA0311', 1, 2, 0, 6, 0 )
-
-    RETURN
-
-    ! Now reverse order of factors
-40  CONTINUE
+    k    = 0
+    l    = 1
+
+    DO
+
+       IF ( MOD( nu, ifac ) == 0 )  THEN
+          k = k + 1
+          jfax(k) = ifac
+          IF ( ifac == 8  .AND.  k /= 1 )  THEN
+             jfax(1) = 8
+             jfax(k) = 6
+          ENDIF
+          nu = nu/ifac
+          IF ( nu == 1   )  EXIT
+          IF ( ifac /= 8 )  CYCLE
+       ENDIF
+
+       l = l + 1
+       ifac = lfax(l)
+       IF (ifac < 2 )  THEN
+          message_string = 'number of gridpoints along x or/and y ' //                             &
+                           'contain illegal  factors' //                                           &
+                           '&only factors 2, 3, 5 are allowed'
+          CALL message( 'temperton_fft', 'PA0311', 1, 2, 0, 6, 0 )
+          RETURN
+       ENDIF
+
+    ENDDO
+!
+!-- Now reverse order of factors
     nfax = k
     ifax(1) = nfax
-    DO i = 1, nfax
+    DO  i = 1, nfax
        ifax(nfax+2-i) = jfax(i)
-    END DO
+    ENDDO
     ifax(10) = n
-    RETURN
-  END SUBROUTINE set99
+
+ END SUBROUTINE set99
 
  END MODULE temperton_fft
Index: palm/trunk/SOURCE/wind_turbine_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/wind_turbine_model_mod.f90	(revision 3724)
+++ palm/trunk/SOURCE/wind_turbine_model_mod.f90	(revision 3725)
@@ -26,4 +26,7 @@
 ! -----------------
 ! $Id$
+! unused variables removed
+! 
+! 3685 2019-01-21 01:02:11Z knoop
 ! Some interface calls moved to module_interface + cleanup
 ! 
@@ -237,5 +240,5 @@
     REAL(wp), DIMENSION(1:100) ::  rnac             = 0.0_wp  !< nacelle diameter [m]
     REAL(wp), DIMENSION(1:100) ::  rr              = 63.0_wp  !< rotor radius [m]
-    REAL(wp), DIMENSION(1:100) ::  turb_cd_nacelle = 0.85_wp  !< drag coefficient for nacelle
+!    REAL(wp), DIMENSION(1:100) ::  turb_cd_nacelle = 0.85_wp  !< drag coefficient for nacelle
     REAL(wp), DIMENSION(1:100) ::  turb_cd_tower    = 1.2_wp  !< drag coefficient for tower
 
@@ -306,5 +309,4 @@
     REAL(wp) ::  eps_min          !< 
     REAL(wp) ::  eps_min2         !< 
-    REAL(wp) ::  sqrt_arg         !< 
 
 !
@@ -519,6 +521,4 @@
        IMPLICIT NONE
        
-       INTEGER(iwp) ::  ierrn       !<
-
        CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file 
 
@@ -531,6 +531,7 @@
                                   rnac, rr, segment_length, segment_width,     &
                                   slope2, speed_control, tilt, time_turbine_on,&
-                                  turb_cd_nacelle, turb_cd_tower, pitch_rate,  &
+                                  turb_cd_tower, pitch_rate,                   &
                                   yaw_control, yaw_speed, tl_cor
+!                                  , turb_cd_nacelle
                                   
        NAMELIST /wind_turbine_parameters/                                      &
@@ -543,6 +544,7 @@
                                   rnac, rr, segment_length, segment_width,     &
                                   slope2, speed_control, tilt, time_turbine_on,&
-                                  turb_cd_nacelle, turb_cd_tower, pitch_rate,  &
+                                  turb_cd_tower, pitch_rate,                   &
                                   yaw_control, yaw_speed, tl_cor
+!                                  , turb_cd_nacelle
 !
 !--    Try to find wind turbine model package
@@ -980,12 +982,4 @@
        INTEGER(iwp) ::  tower_n      !<
        INTEGER(iwp) ::  tower_s      !<
-!
-!--    Help variables for the calulaction of the nacelle drag
-       INTEGER(iwp) ::  i_ip         !<
-       INTEGER(iwp) ::  i_ipg        !<
-       
-       REAL(wp) ::  yvalue               
-       REAL(wp) ::  dy_int           !< 
-       REAL(wp) ::  dz_int           !<
        
        REAL(wp), DIMENSION(:,:), ALLOCATABLE :: circle_points  !<
