Changeset 3725


Ignore:
Timestamp:
Feb 7, 2019 10:11:02 AM (3 years ago)
Author:
raasch
Message:

modifications to avoid compiler warnings about unused variables, temperton-fft: GOTO statements replaced, file re-formatted corresponding to coding standards, ssh-calls for compilations on remote systems modified to avoid output of login messages on specific systems changed again (palmbuild, reverted as before r3549), error messages for failed restarts extended (palmrun)

Location:
palm/trunk
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SCRIPTS/palmbuild

    r3549 r3725  
    2727# -----------------
    2828# $Id$
     29# ssh-calls for compilations on remote systems modified to avoid output
     30# of login messages on specific systems changed again (reverted as before r3549)
     31#
     32# 3549 2018-11-21 15:44:44Z raasch
    2933# ssh-calls for compilations on remote systems modified to avoid output
    3034# of login messages on specific systems
     
    571575       fi
    572576       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\" "
    573        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
    574 ###       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
     577###       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
     578       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
    575579
    576580       if [[ $(grep -c MAKE_ERROR ${configuration_identifier}_last_make_protocol) != 0 ]]
     
    611615       fi
    612616       make_call_string="make  $make_options  PROG=$program_name  F90=$compiler_name  COPT=\"$cpp_options\"  F90FLAGS=\"$compiler_options\"  LDFLAGS=\"$linker_options\" "
    613        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
    614 ###       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
     617###       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
     618       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
    615619
    616620       if [[ $(grep -c MAKE_ERROR ${configuration_identifier}_last_make_protocol) != 0 ]]
  • palm/trunk/SCRIPTS/palmrun

    r3682 r3725  
    2727# -----------------
    2828# $Id$
     29# error messages for failed restarts extended
     30#
     31# 3682 2019-01-18 17:01:54Z knoop
    2932# ssh-call for submitting batch jobs on remote systems modified again to avoid output
    3033# of login messages on specific systems
     
    27992802
    28002803          # CHECK, IF RESTART JOB HAS BEEN STARTED
    2801        if [[ $( grep -c "+++ palmrun killed"  palmrun_restart.log ) != 0 ]]
     2804       if [[ $( grep -c "+++ palmrun crashed"  palmrun_restart.log ) != 0 ]]
    28022805       then
    28032806          printf "\n$dashes\n  +++ creating restart run failed \n"
     
    28052808          rm  palmrun_restart.log
    28062809          exit
     2810       elif [[ $( grep -c "*** palmrun finished"  palmrun_restart.log ) != 1 ]]
     2811       then
     2812          printf "\n$dashes\n  +++ creating restart run failed, probably due to network problems\n"
     2813          locat=create_restart
     2814          rm  palmrun_restart.log
     2815          exit       
    28072816       else
    28082817          printf "\n$dashes\n  *** restart run initiated \n"
  • palm/trunk/SOURCE/buoyancy.f90

    r3655 r3725  
    2525! -----------------
    2626! $Id$
     27! unused variables removed
     28!
     29! 3655 2019-01-07 16:51:22Z knoop
    2730! nopointer option removed
    2831!
     
    160163
    161164       USE indices,                                                            &
    162            ONLY:  nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nzb,       &
    163                   nzt
     165           ONLY:  nxl, nxlu, nxr, nyn, nys, nzb, nzt
    164166
    165167
     
    251253
    252254       USE indices,                                                            &
    253            ONLY:  nxlg, nxrg, nyng, nysg, nzb, nzt
     255           ONLY:  nzb, nzt
    254256
    255257       IMPLICIT NONE
  • palm/trunk/SOURCE/data_log.f90

    r3655 r3725  
    2525! -----------------
    2626! $Id$
     27! preprocessor directives removed to avoid compiler warnings
     28!
     29! 3655 2019-01-07 16:51:22Z knoop
    2730! Corrected "Former revisions" section
    2831!
     
    6164 SUBROUTINE data_log( array, i1, i2, j1, j2, k1, k2 )
    6265 
    63 #if defined( __logging )
    64 
    6566    USE control_parameters,                                                    &
    6667        ONLY:  log_message, simulated_time
     
    9899    WRITE ( 20 )  array
    99100
    100 #endif
    101101 END SUBROUTINE data_log
    102102
     
    110110 
    111111 SUBROUTINE data_log_2d( array, i1, i2, j1, j2)
    112 
    113 #if defined( __logging )
    114112
    115113    USE control_parameters,                                                    &
     
    146144    WRITE ( 20 )  array
    147145
    148 #endif
    149146 END SUBROUTINE data_log_2d
    150147
     
    158155 
    159156 SUBROUTINE data_log_2d_int( array, i1, i2, j1, j2)
    160 
    161 #if defined( __logging )
    162157
    163158    USE control_parameters,                                                    &
     
    194189    WRITE ( 20 )  array
    195190
    196 #endif
    197191 END SUBROUTINE data_log_2d_int
  • palm/trunk/SOURCE/gust_mod.f90

    r3685 r3725  
    2525! -----------------
    2626! $Id$
     27! dummy statement modified to avoid compiler warnings about unused variables
     28!
     29! 3685 2019-01-21 01:02:11Z knoop
    2730! Some interface calls moved to module_interface + cleanup
    2831!
     
    445448!
    446449!--    Next line is just to avoid compiler warnings about unused variables. You may remove it.
    447        IF ( found .AND. two_d )  idum = av + LEN( variable ) + LEN( grid ) + local_pf(nxl,nys,nzb_do) + fill_value
     450       IF ( found .AND. two_d )  THEN
     451          idum = av + LEN( variable ) + LEN( grid // mode ) + local_pf(nxl,nys,nzb_do) + fill_value
     452       ENDIF
    448453
    449454    END SUBROUTINE gust_data_output_2d
  • palm/trunk/SOURCE/poismg_mod.f90

    r3655 r3725  
    2525! -----------------
    2626! $Id$
     27! unused subroutine removed
     28!
     29! 3655 2019-01-07 16:51:22Z knoop
    2730! unnecessary check eliminated
    2831!
     
    143146    INTERFACE sort_k_to_even_odd_blocks
    144147       MODULE PROCEDURE sort_k_to_even_odd_blocks
    145        MODULE PROCEDURE sort_k_to_even_odd_blocks_int
     148!       MODULE PROCEDURE sort_k_to_even_odd_blocks_int
    146149       MODULE PROCEDURE sort_k_to_even_odd_blocks_1d
    147150    END INTERFACE sort_k_to_even_odd_blocks
     
    11521155!> Version for 2D-INTEGER arrays
    11531156!------------------------------------------------------------------------------!
    1154     SUBROUTINE sort_k_to_even_odd_blocks_int( i_mg , glevel )
    1155 
    1156 
    1157        USE indices,                                                            &
    1158            ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
    1159 
    1160        IMPLICIT NONE
    1161 
    1162        INTEGER(iwp), INTENT(IN) ::  glevel  !< grid level
    1163 
    1164        INTEGER(iwp), DIMENSION(nzb:nzt_mg(glevel)+1,                           &
    1165                                nys_mg(glevel)-1:nyn_mg(glevel)+1,              &
    1166                                nxl_mg(glevel)-1:nxr_mg(glevel)+1) ::           &
    1167                                     i_mg    !< array to be sorted
     1157!    SUBROUTINE sort_k_to_even_odd_blocks_int( i_mg , glevel )
     1158!
     1159!
     1160!       USE indices,                                                            &
     1161!           ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
     1162!
     1163!       IMPLICIT NONE
     1164!
     1165!       INTEGER(iwp), INTENT(IN) ::  glevel  !< grid level
     1166!
     1167!       INTEGER(iwp), DIMENSION(nzb:nzt_mg(glevel)+1,                           &
     1168!                               nys_mg(glevel)-1:nyn_mg(glevel)+1,              &
     1169!                               nxl_mg(glevel)-1:nxr_mg(glevel)+1) ::           &
     1170!                                    i_mg    !< array to be sorted
    11681171!
    11691172!--    Local variables
    1170        INTEGER(iwp) :: i        !< index variabel along x
    1171        INTEGER(iwp) :: j        !< index variable along y
    1172        INTEGER(iwp) :: k        !< index variable along z
    1173        INTEGER(iwp) :: l        !< grid level
    1174        INTEGER(iwp) :: ind      !< index variable along z
    1175        INTEGER(iwp),DIMENSION(nzb:nzt_mg(glevel)+1) ::  tmp  !< temporary odd-even sorted array
    1176 
    1177 
    1178        CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'start' )
    1179 
    1180        l = glevel
    1181        ind_even_odd = even_odd_level(l)
    1182 
    1183        DO  i = nxl_mg(l)-1, nxr_mg(l)+1
    1184           DO  j = nys_mg(l)-1, nyn_mg(l)+1
    1185 
     1173!       INTEGER(iwp) :: i        !< index variabel along x
     1174!       INTEGER(iwp) :: j        !< index variable along y
     1175!       INTEGER(iwp) :: k        !< index variable along z
     1176!       INTEGER(iwp) :: l        !< grid level
     1177!       INTEGER(iwp) :: ind      !< index variable along z
     1178!       INTEGER(iwp),DIMENSION(nzb:nzt_mg(glevel)+1) ::  tmp  !< temporary odd-even sorted array
     1179!
     1180!
     1181!       CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'start' )
     1182!
     1183!       l = glevel
     1184!       ind_even_odd = even_odd_level(l)
     1185!
     1186!       DO  i = nxl_mg(l)-1, nxr_mg(l)+1
     1187!          DO  j = nys_mg(l)-1, nyn_mg(l)+1
     1188!
    11861189!
    11871190!--          Sort the data with even k index
    1188              ind = nzb-1
    1189              DO  k = nzb, nzt_mg(l), 2
    1190                 ind = ind + 1
    1191                 tmp(ind) = i_mg(k,j,i)
    1192              ENDDO
     1191!             ind = nzb-1
     1192!             DO  k = nzb, nzt_mg(l), 2
     1193!                ind = ind + 1
     1194!                tmp(ind) = i_mg(k,j,i)
     1195!             ENDDO
    11931196!
    11941197!--          Sort the data with odd k index
    1195              DO  k = nzb+1, nzt_mg(l)+1, 2
    1196                 ind = ind + 1
    1197                 tmp(ind) = i_mg(k,j,i)
    1198              ENDDO
    1199 
    1200              i_mg(:,j,i) = tmp
    1201 
    1202           ENDDO
    1203        ENDDO
    1204 
    1205        CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'stop' )
    1206 
    1207     END SUBROUTINE sort_k_to_even_odd_blocks_int
     1198!             DO  k = nzb+1, nzt_mg(l)+1, 2
     1199!                ind = ind + 1
     1200!                tmp(ind) = i_mg(k,j,i)
     1201!             ENDDO
     1202!
     1203!             i_mg(:,j,i) = tmp
     1204!
     1205!          ENDDO
     1206!       ENDDO
     1207!
     1208!       CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'stop' )
     1209!
     1210!    END SUBROUTINE sort_k_to_even_odd_blocks_int
    12081211
    12091212
  • palm/trunk/SOURCE/temperton_fft_mod.f90

    r3241 r3725  
    99! -----------------
    1010! $Id$
     11! GOTO statements replaced, file re-formatted corresponding to coding standards
     12!
     13! 3241 2018-09-12 15:02:00Z raasch
    1114! unused variables removed
    1215!
     
    2528! Module renamed
    2629!
    27 !
    2830! 1682 2015-10-07 23:56:08Z knoop
    2931! Code annotations made doxygen readable
     
    7072! Description:
    7173! ------------
    72 !> Calls fortran-versions of fft's.
     74!> Calls Fortran-versions of fft's.
    7375!>
    7476!> Method:
     
    115117!> dimension a(n),work(n),trigs(n),ifax(1)
    116118!------------------------------------------------------------------------------!
    117   SUBROUTINE fft991cy(a,work,trigs,ifax,inc,jump,n,lot,isign)
     119 SUBROUTINE fft991cy( a, work, trigs, ifax, inc, jump, n, lot, isign )
    118120
    119121    USE kinds
     
    121123    IMPLICIT NONE
    122124
    123     !  Scalar arguments
    124125    INTEGER(iwp) ::  inc   !<
    125126    INTEGER(iwp) ::  isign !<
     
    128129    INTEGER(iwp) ::  n     !<
    129130
    130     !  Array arguments
    131131    REAL(wp)     ::  a(*)     !<
    132132    REAL(wp)     ::  trigs(*) !<
     
    134134    INTEGER(iwp) ::  ifax(*)  !<
    135135
    136     !  Local scalars:
    137136    INTEGER(iwp) ::  i      !<
    138137    INTEGER(iwp) ::  ia     !<
     
    156155    INTEGER(iwp) ::  nx     !<
    157156
    158     !  Intrinsic functions
    159 !    INTRINSIC MOD
    160 
    161 
    162     !  Executable statements
    163 
    164     IF (ifax(10)/=n) CALL set99(trigs,ifax,n)
     157
     158    IF ( ifax(10) /= n )  CALL set99( trigs, ifax, n )
    165159    nfax = ifax(1)
    166     nx = n + 1
    167     IF (MOD(n,2)==1) nx = n
    168     nblox = 1 + (lot-1)/nfft
    169     nvex = lot - (nblox-1)*nfft
    170     IF (isign==-1) GO TO 50
    171 
    172     ! isign=+1, spectral to gridpoint transform
    173 
    174     istart = 1
    175     DO nb = 1, nblox
    176        ia = istart
    177        i = istart
    178 !DIR$ IVDEP
    179 !CDIR NODEP
    180 !OCL NOVREC
    181        DO j = 1, nvex
    182           a(i+inc) = 0.5_wp*a(i)
    183           i = i + jump
    184        END DO
    185        IF (MOD(n,2)==1) GO TO 10
    186        i = istart + n*inc
    187        DO j = 1, nvex
    188           a(i) = 0.5_wp*a(i)
    189           i = i + jump
    190        END DO
    191 10     CONTINUE
    192        ia = istart + inc
    193        la = 1
    194        igo = + 1
    195 
    196        DO k = 1, nfax
    197           ifac = ifax(k+1)
    198           ierr = -1
    199           IF (igo==-1) GO TO 20
    200           CALL rpassm(a(ia),a(ia+la*inc),work(1),work(ifac*la+1),trigs,inc,1, &
    201                &          jump,nx,nvex,n,ifac,la,ierr)
    202           GO TO 30
    203 20        CONTINUE
    204           CALL rpassm(work(1),work(la+1),a(ia),a(ia+ifac*la*inc),trigs,1,inc,nx, &
    205                &          jump,nvex,n,ifac,la,ierr)
    206 30        CONTINUE
    207           IF (ierr/=0) GO TO 100
    208           la = ifac*la
    209           igo = -igo
     160    nx   = n + 1
     161    IF ( MOD(n,2) == 1 )  nx = n
     162    nblox = 1 + (lot-1) / nfft
     163    nvex = lot - (nblox-1) * nfft
     164
     165
     166    IF ( isign == 1 )  THEN
     167!
     168!--    Backward fft: spectral to gridpoint transform
     169
     170       istart = 1
     171       DO  nb = 1, nblox
     172
    210173          ia = istart
    211        END DO
    212 
    213        ! If necessary, copy results back to a
    214 
    215        IF (MOD(nfax,2)==0) GO TO 40
    216        ibase = 1
    217        jbase = ia
    218        DO jj = 1, nvex
    219           i = ibase
    220           j = jbase
    221           DO ii = 1, n
    222              a(j) = work(i)
    223              i = i + 1
    224              j = j + inc
    225           END DO
    226           ibase = ibase + nx
    227           jbase = jbase + jump
    228        END DO
    229 40     CONTINUE
    230 
    231        ! Fill in zeros at end
    232 
    233        ix = istart + n*inc
    234 !DIR$ IVDEP
    235 !CDIR NODEP
    236 !OCL NOVREC
    237        DO j = 1, nvex
    238           a(ix) = 0.0_wp
    239           a(ix+inc) = 0.0_wp
    240           ix = ix + jump
    241        END DO
    242 
    243        istart = istart + nvex*jump
    244        nvex = nfft
    245     END DO
    246     RETURN
    247 
    248     ! isign=-1, gridpoint to spectral transform
    249 
    250 50  CONTINUE
    251     istart = 1
    252     DO nb = 1, nblox
    253        ia = istart
    254        la = n
    255        igo = + 1
    256 
    257        DO k = 1, nfax
    258           ifac = ifax(nfax+2-k)
    259           la = la/ifac
    260           ierr = -1
    261           IF (igo==-1) GO TO 60
    262           CALL qpassm(a(ia),a(ia+ifac*la*inc),work(1),work(la+1),trigs,inc,1, &
    263                &          jump,nx,nvex,n,ifac,la,ierr)
    264           GO TO 70
    265 60        CONTINUE
    266           CALL qpassm(work(1),work(ifac*la+1),a(ia),a(ia+la*inc),trigs,1,inc,nx, &
    267                &          jump,nvex,n,ifac,la,ierr)
    268 70        CONTINUE
    269           IF (ierr/=0) GO TO 100
    270           igo = -igo
     174          i = istart
     175          !DIR$ IVDEP
     176          !CDIR NODEP
     177          !OCL NOVREC
     178          DO  j = 1, nvex
     179             a(i+inc) = 0.5_wp * a(i)
     180             i = i + jump
     181          ENDDO
     182          IF ( MOD(n,2) /= 1 )  THEN
     183             i = istart + n * inc
     184             DO  j = 1, nvex
     185                a(i) = 0.5_wp * a(i)
     186                i = i + jump
     187             ENDDO
     188          ENDIF
    271189          ia = istart + inc
    272        END DO
    273 
    274        ! If necessary, copy results back to a
    275 
    276        IF (MOD(nfax,2)==0) GO TO 80
    277        ibase = 1
    278        jbase = ia
    279        DO jj = 1, nvex
    280           i = ibase
    281           j = jbase
    282           DO ii = 1, n
    283              a(j) = work(i)
    284              i = i + 1
    285              j = j + inc
    286           END DO
    287           ibase = ibase + nx
    288           jbase = jbase + jump
    289        END DO
    290 80     CONTINUE
    291 
    292        ! Shift a(0) & fill in zero imag parts
    293 
    294        ix = istart
    295 !DIR$ IVDEP
    296 !CDIR NODEP
    297 !OCL NOVREC
    298        DO j = 1, nvex
    299           a(ix) = a(ix+inc)
    300           a(ix+inc) = 0.0_wp
    301           ix = ix + jump
    302        END DO
    303        IF (MOD(n,2)==1) GO TO 90
    304        iz = istart + (n+1)*inc
    305        DO j = 1, nvex
    306           a(iz) = 0.0_wp
    307           iz = iz + jump
    308        END DO
    309 90     CONTINUE
    310 
    311        istart = istart + nvex*jump
    312        nvex = nfft
    313     END DO
    314     RETURN
    315 
    316     ! Error messages
    317 
    318 100 CONTINUE
    319 
    320     SELECT CASE (ierr)
    321     CASE (:-1)
    322        WRITE (nout,'(A,I5,A)') ' Vector length =',nvex,', greater than nfft'
    323     CASE (0)
    324        WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', not catered for'
    325     CASE (1:)
    326        WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', only catered for if la*ifac=n'
    327     END SELECT
    328 
    329     RETURN
    330   END SUBROUTINE fft991cy
     190          la = 1
     191          igo = + 1
     192
     193          DO  k = 1, nfax
     194
     195             ifac = ifax(k+1)
     196             ierr = -1
     197             IF ( igo /= -1 )  THEN
     198                CALL rpassm( a(ia), a(ia+la*inc), work(1), work(ifac*la+1), trigs, inc, 1, jump,   &
     199                             nx, nvex, n, ifac, la, ierr )
     200             ELSE
     201                CALL rpassm( work(1), work(la+1), a(ia), a(ia+ifac*la*inc), trigs, 1, inc, nx,     &
     202                             jump, nvex, n, ifac, la, ierr )
     203             ENDIF
     204!
     205!--          Following messages shouldn't appear in PALM applications
     206             IF ( ierr /= 0 )  THEN
     207                SELECT CASE (ierr)
     208                   CASE (:-1)
     209                      WRITE (nout,'(A,I5,A)') ' Vector length =',nvex,', greater than nfft'
     210                   CASE (0)
     211                      WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', not catered for'
     212                   CASE (1:)
     213                      WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', only catered for if la*ifac=n'
     214                END SELECT
     215                RETURN
     216             ENDIF
     217
     218             la  = ifac * la
     219             igo = -igo
     220             ia  = istart
     221
     222          ENDDO
     223!
     224!--       If necessary, copy results back to a
     225          IF ( MOD(nfax,2) /= 0 )  THEN
     226             ibase = 1
     227             jbase = ia
     228             DO  jj = 1, nvex
     229                i = ibase
     230                j = jbase
     231                DO  ii = 1, n
     232                   a(j) = work(i)
     233                   i = i + 1
     234                   j = j + inc
     235                ENDDO
     236                ibase = ibase + nx
     237                jbase = jbase + jump
     238             ENDDO
     239          ENDIF
     240!
     241!--       Fill in zeros at end
     242          ix = istart + n*inc
     243          !DIR$ IVDEP
     244          !CDIR NODEP
     245          !OCL NOVREC
     246          DO  j = 1, nvex
     247             a(ix) = 0.0_wp
     248             a(ix+inc) = 0.0_wp
     249             ix = ix + jump
     250          ENDDO
     251          istart = istart + nvex*jump
     252          nvex = nfft
     253
     254       ENDDO
     255
     256    ELSEIF ( isign == -1 )  THEN
     257!
     258!--    Forward fft: gridpoint to spectral transform
     259       istart = 1
     260       DO  nb = 1, nblox
     261          ia  = istart
     262          la  = n
     263          igo = + 1
     264
     265          DO  k = 1, nfax
     266
     267             ifac = ifax(nfax+2-k)
     268             la = la / ifac
     269             ierr = -1
     270             IF ( igo /= -1 )  THEN
     271                CALL qpassm( a(ia), a(ia+ifac*la*inc), work(1), work(la+1), trigs, inc, 1, jump,   &
     272                             nx, nvex, n, ifac, la, ierr )
     273             ELSE
     274                CALL qpassm( work(1), work(ifac*la+1), a(ia), a(ia+la*inc), trigs, 1, inc, nx,     &
     275                             jump, nvex, n, ifac, la, ierr )
     276             ENDIF
     277!
     278!--          Following messages shouldn't appear in PALM applications
     279             IF ( ierr /= 0 )  THEN
     280                SELECT CASE (ierr)
     281                   CASE (:-1)
     282                      WRITE (nout,'(A,I5,A)') ' Vector length =',nvex,', greater than nfft'
     283                   CASE (0)
     284                      WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', not catered for'
     285                   CASE (1:)
     286                      WRITE (nout,'(A,I3,A)') ' Factor =',ifac,', only catered for if la*ifac=n'
     287                END SELECT
     288                RETURN
     289             ENDIF
     290
     291             igo = -igo
     292             ia = istart + inc
     293
     294          ENDDO
     295!
     296!--       If necessary, copy results back to a
     297          IF ( MOD(nfax,2) /= 0 )  THEN
     298             ibase = 1
     299             jbase = ia
     300             DO  jj = 1, nvex
     301                i = ibase
     302                j = jbase
     303                DO  ii = 1, n
     304                   a(j) = work(i)
     305                   i = i + 1
     306                   j = j + inc
     307                ENDDO
     308                ibase = ibase + nx
     309                jbase = jbase + jump
     310             ENDDO
     311          ENDIF
     312!
     313!--       Shift a(0) and fill in zero imag parts
     314          ix = istart
     315          !DIR$ IVDEP
     316          !CDIR NODEP
     317          !OCL NOVREC
     318          DO  j = 1, nvex
     319             a(ix) = a(ix+inc)
     320             a(ix+inc) = 0.0_wp
     321             ix = ix + jump
     322          ENDDO
     323          IF ( MOD(n,2) /= 1 )  THEN
     324             iz = istart + (n+1) * inc
     325             DO  j = 1, nvex
     326                a(iz) = 0.0_wp
     327                iz = iz + jump
     328             ENDDO
     329          ENDIF
     330          istart = istart + nvex * jump
     331          nvex = nfft
     332
     333       ENDDO
     334
     335    ENDIF
     336
     337 END SUBROUTINE fft991cy
    331338
    332339!------------------------------------------------------------------------------!
     
    357364!>         3 - ifac only catered for if la=n/ifac
    358365!------------------------------------------------------------------------------!
    359   SUBROUTINE qpassm(a,b,c,d,trigs,inc1,inc2,inc3,inc4,lot,n,ifac,la,ierr)
     366 SUBROUTINE qpassm( a, b, c, d, trigs, inc1, inc2, inc3, inc4, lot, n, ifac, la, ierr )
    360367
    361368    USE kinds
     
    363370    IMPLICIT NONE
    364371
    365     !  Scalar arguments
    366372    INTEGER(iwp) ::  ierr !<
    367373    INTEGER(iwp) ::  ifac !<
     
    374380    INTEGER(iwp) ::  n    !<
    375381
    376     !  Array arguments
    377     ! REAL :: a(n),b(n),c(n),d(n),trigs(n)
     382!
     383!-- Arrays are dimensioned with n
    378384    REAL(wp) ::  a(*)     !<
    379385    REAL(wp) ::  b(*)     !<
     
    382388    REAL(wp) ::  trigs(*) !<
    383389 
    384     !  Local scalars:
    385390    REAL(wp) ::  a0     !<
    386391    REAL(wp) ::  a1     !<
     
    430435    INTEGER(iwp) ::  ia    !<
    431436    INTEGER(iwp) ::  ib    !<
    432     INTEGER(iwp) ::  ibad  !<
    433437    INTEGER(iwp) ::  ibase !<
    434438    INTEGER(iwp) ::  ic    !<
     
    461465    INTEGER(iwp) ::  m     !<
    462466
    463     !  Intrinsic functions
    464 !    INTRINSIC REAL, SQRT
    465 
    466     !  Data statements
    467     DATA sin36/0.587785252292473_wp/, sin72/0.951056516295154_wp/, &
    468          &      qrt5/0.559016994374947_wp/, sin60/0.866025403784437_wp/
    469 
    470 
    471     !  Executable statements
    472 
    473     m = n/ifac
    474     iink = la*inc1
    475     jink = la*inc2
    476     ijump = (ifac-1)*iink
    477     kstop = (n-ifac)/(2*ifac)
    478 
    479     ibad = 1
    480     IF (lot>nfft) GO TO 180
     467
     468    DATA  sin36/0.587785252292473_wp/, sin72/0.951056516295154_wp/,                                &
     469          qrt5/0.559016994374947_wp/,  sin60/0.866025403784437_wp/
     470
     471
     472    ierr = 0
     473
     474    IF ( lot > nfft )  THEN
     475       ierr = 1
     476       RETURN
     477    ENDIF
     478
     479    m = n / ifac
     480    iink  = la * inc1
     481    jink  = la * inc2
     482    ijump = (ifac-1) * iink
     483    kstop = ( n-ifac ) / ( 2*ifac )
     484
    481485    ibase = 0
    482486    jbase = 0
    483487    igo = ifac - 1
    484     IF (igo==7) igo = 6
    485     ibad = 2
    486     IF (igo<1 .OR. igo>6) GO TO 180
    487     GO TO (10,40,70,100,130,160) igo
    488 
    489     ! Coding for factor 2
    490 
    491 10  CONTINUE
    492     ia = 1
    493     ib = ia + iink
    494     ja = 1
    495     jb = ja + (2*m-la)*inc2
    496 
    497     IF (la==m) GO TO 30
    498 
    499     DO l = 1, la
    500        i = ibase
    501        j = jbase
    502 !DIR$ IVDEP
    503 !CDIR NODEP
    504 !OCL NOVREC
    505        DO ijk = 1, lot
    506           c(ja+j) = a(ia+i) + a(ib+i)
    507           c(jb+j) = a(ia+i) - a(ib+i)
    508           i = i + inc3
    509           j = j + inc4
    510        END DO
    511        ibase = ibase + inc1
    512        jbase = jbase + inc2
    513     END DO
    514     ja = ja + jink
    515     jink = 2*jink
    516     jb = jb - jink
    517     ibase = ibase + ijump
    518     ijump = 2*ijump + iink
    519     IF (ja==jb) GO TO 20
    520     DO k = la, kstop, la
    521        kb = k + k
    522        c1 = trigs(kb+1)
    523        s1 = trigs(kb+2)
    524        jbase = 0
    525        DO l = 1, la
    526           i = ibase
    527           j = jbase
    528 !DIR$ IVDEP
    529 !CDIR NODEP
    530 !OCL NOVREC
    531           DO ijk = 1, lot
    532              c(ja+j) = a(ia+i) + (c1*a(ib+i)+s1*b(ib+i))
    533              c(jb+j) = a(ia+i) - (c1*a(ib+i)+s1*b(ib+i))
    534              d(ja+j) = (c1*b(ib+i)-s1*a(ib+i)) + b(ia+i)
    535              d(jb+j) = (c1*b(ib+i)-s1*a(ib+i)) - b(ia+i)
    536              i = i + inc3
    537              j = j + inc4
    538           END DO
    539           ibase = ibase + inc1
    540           jbase = jbase + inc2
    541        END DO
    542        ibase = ibase + ijump
    543        ja = ja + jink
    544        jb = jb - jink
    545     END DO
    546     IF (ja>jb) GO TO 170
    547 20  CONTINUE
    548     jbase = 0
    549     DO l = 1, la
    550        i = ibase
    551        j = jbase
    552 !DIR$ IVDEP
    553 !CDIR NODEP
    554 !OCL NOVREC
    555        DO ijk = 1, lot
    556           c(ja+j) = a(ia+i)
    557           d(ja+j) = -a(ib+i)
    558           i = i + inc3
    559           j = j + inc4
    560        END DO
    561        ibase = ibase + inc1
    562        jbase = jbase + inc2
    563     END DO
    564 
    565     GO TO 170
    566 30  CONTINUE
    567     z = 1.0_wp/REAL(n,KIND=wp)
    568     DO l = 1, la
    569        i = ibase
    570        j = jbase
    571 !DIR$ IVDEP
    572 !CDIR NODEP
    573 !OCL NOVREC
    574        DO ijk = 1, lot
    575           c(ja+j) = z*(a(ia+i)+a(ib+i))
    576           c(jb+j) = z*(a(ia+i)-a(ib+i))
    577           i = i + inc3
    578           j = j + inc4
    579        END DO
    580        ibase = ibase + inc1
    581        jbase = jbase + inc2
    582     END DO
    583     GO TO 170
    584 
    585     ! Coding for factor 3
    586 
    587 40  CONTINUE
    588     ia = 1
    589     ib = ia + iink
    590     ic = ib + iink
    591     ja = 1
    592     jb = ja + (2*m-la)*inc2
    593     jc = jb
    594 
    595     IF (la==m) GO TO 60
    596 
    597     DO l = 1, la
    598        i = ibase
    599        j = jbase
    600 !DIR$ IVDEP
    601 !CDIR NODEP
    602 !OCL NOVREC
    603        DO ijk = 1, lot
    604           c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
    605           c(jb+j) = a(ia+i) - 0.5_wp*(a(ib+i)+a(ic+i))
    606           d(jb+j) = sin60*(a(ic+i)-a(ib+i))
    607           i = i + inc3
    608           j = j + inc4
    609        END DO
    610        ibase = ibase + inc1
    611        jbase = jbase + inc2
    612     END DO
    613     ja = ja + jink
    614     jink = 2*jink
    615     jb = jb + jink
    616     jc = jc - jink
    617     ibase = ibase + ijump
    618     ijump = 2*ijump + iink
    619     IF (ja==jc) GO TO 50
    620     DO k = la, kstop, la
    621        kb = k + k
    622        kc = kb + kb
    623        c1 = trigs(kb+1)
    624        s1 = trigs(kb+2)
    625        c2 = trigs(kc+1)
    626        s2 = trigs(kc+2)
    627        jbase = 0
    628        DO l = 1, la
    629           i = ibase
    630           j = jbase
    631 !DIR$ IVDEP
    632 !CDIR NODEP
    633 !OCL NOVREC
    634           DO ijk = 1, lot
    635              a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c2*a(ic+i)+s2*b(ic+i))
    636              b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c2*b(ic+i)-s2*a(ic+i))
    637              a2 = a(ia+i) - 0.5_wp*a1
    638              b2 = b(ia+i) - 0.5_wp*b1
    639              a3 = sin60*((c1*a(ib+i)+s1*b(ib+i))-(c2*a(ic+i)+s2*b(ic+i)))
    640              b3 = sin60*((c1*b(ib+i)-s1*a(ib+i))-(c2*b(ic+i)-s2*a(ic+i)))
    641              c(ja+j) = a(ia+i) + a1
    642              d(ja+j) = b(ia+i) + b1
    643              c(jb+j) = a2 + b3
    644              d(jb+j) = b2 - a3
    645              c(jc+j) = a2 - b3
    646              d(jc+j) = -(b2+a3)
    647              i = i + inc3
    648              j = j + inc4
    649           END DO
    650           ibase = ibase + inc1
    651           jbase = jbase + inc2
    652        END DO
    653        ibase = ibase + ijump
    654        ja = ja + jink
    655        jb = jb + jink
    656        jc = jc - jink
    657     END DO
    658     IF (ja>jc) GO TO 170
    659 50  CONTINUE
    660     jbase = 0
    661     DO l = 1, la
    662        i = ibase
    663        j = jbase
    664 !DIR$ IVDEP
    665 !CDIR NODEP
    666 !OCL NOVREC
    667        DO ijk = 1, lot
    668           c(ja+j) = a(ia+i) + 0.5_wp*(a(ib+i)-a(ic+i))
    669           d(ja+j) = -sin60*(a(ib+i)+a(ic+i))
    670           c(jb+j) = a(ia+i) - (a(ib+i)-a(ic+i))
    671           i = i + inc3
    672           j = j + inc4
    673        END DO
    674        ibase = ibase + inc1
    675        jbase = jbase + inc2
    676     END DO
    677 
    678     GO TO 170
    679 60  CONTINUE
    680     z = 1.0_wp/REAL(n,KIND=wp)
    681     zsin60 = z*sin60
    682     DO l = 1, la
    683        i = ibase
    684        j = jbase
    685 !DIR$ IVDEP
    686 !CDIR NODEP
    687 !OCL NOVREC
    688        DO ijk = 1, lot
    689           c(ja+j) = z*(a(ia+i)+(a(ib+i)+a(ic+i)))
    690           c(jb+j) = z*(a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))
    691           d(jb+j) = zsin60*(a(ic+i)-a(ib+i))
    692           i = i + inc3
    693           j = j + inc4
    694        END DO
    695        ibase = ibase + inc1
    696        jbase = jbase + inc2
    697     END DO
    698     GO TO 170
    699 
    700     ! Coding for factor 4
    701 
    702 70  CONTINUE
    703     ia = 1
    704     ib = ia + iink
    705     ic = ib + iink
    706     id = ic + iink
    707     ja = 1
    708     jb = ja + (2*m-la)*inc2
    709     jc = jb + 2*m*inc2
    710     jd = jb
    711 
    712     IF (la==m) GO TO 90
    713 
    714     DO l = 1, la
    715        i = ibase
    716        j = jbase
    717 !DIR$ IVDEP
    718 !CDIR NODEP
    719 !OCL NOVREC
    720        DO ijk = 1, lot
    721           c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))
    722           c(jc+j) = (a(ia+i)+a(ic+i)) - (a(ib+i)+a(id+i))
    723           c(jb+j) = a(ia+i) - a(ic+i)
    724           d(jb+j) = a(id+i) - a(ib+i)
    725           i = i + inc3
    726           j = j + inc4
    727        END DO
    728        ibase = ibase + inc1
    729        jbase = jbase + inc2
    730     END DO
    731     ja = ja + jink
    732     jink = 2*jink
    733     jb = jb + jink
    734     jc = jc - jink
    735     jd = jd - jink
    736     ibase = ibase + ijump
    737     ijump = 2*ijump + iink
    738     IF (jb==jc) GO TO 80
    739     DO k = la, kstop, la
    740        kb = k + k
    741        kc = kb + kb
    742        kd = kc + kb
    743        c1 = trigs(kb+1)
    744        s1 = trigs(kb+2)
    745        c2 = trigs(kc+1)
    746        s2 = trigs(kc+2)
    747        c3 = trigs(kd+1)
    748        s3 = trigs(kd+2)
    749        jbase = 0
    750        DO l = 1, la
    751           i = ibase
    752           j = jbase
    753 !DIR$ IVDEP
    754 !CDIR NODEP
    755 !OCL NOVREC
    756           DO ijk = 1, lot
    757              a0 = a(ia+i) + (c2*a(ic+i)+s2*b(ic+i))
    758              a2 = a(ia+i) - (c2*a(ic+i)+s2*b(ic+i))
    759              a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c3*a(id+i)+s3*b(id+i))
    760              a3 = (c1*a(ib+i)+s1*b(ib+i)) - (c3*a(id+i)+s3*b(id+i))
    761              b0 = b(ia+i) + (c2*b(ic+i)-s2*a(ic+i))
    762              b2 = b(ia+i) - (c2*b(ic+i)-s2*a(ic+i))
    763              b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c3*b(id+i)-s3*a(id+i))
    764              b3 = (c1*b(ib+i)-s1*a(ib+i)) - (c3*b(id+i)-s3*a(id+i))
    765              c(ja+j) = a0 + a1
    766              c(jc+j) = a0 - a1
    767              d(ja+j) = b0 + b1
    768              d(jc+j) = b1 - b0
    769              c(jb+j) = a2 + b3
    770              c(jd+j) = a2 - b3
    771              d(jb+j) = b2 - a3
    772              d(jd+j) = -(b2+a3)
    773              i = i + inc3
    774              j = j + inc4
    775           END DO
    776           ibase = ibase + inc1
    777           jbase = jbase + inc2
    778        END DO
    779        ibase = ibase + ijump
    780        ja = ja + jink
    781        jb = jb + jink
    782        jc = jc - jink
    783        jd = jd - jink
    784     END DO
    785     IF (jb>jc) GO TO 170
    786 80  CONTINUE
    787     sin45 = SQRT(0.5_wp)
    788     jbase = 0
    789     DO l = 1, la
    790        i = ibase
    791        j = jbase
    792 !DIR$ IVDEP
    793 !CDIR NODEP
    794 !OCL NOVREC
    795        DO ijk = 1, lot
    796           c(ja+j) = a(ia+i) + sin45*(a(ib+i)-a(id+i))
    797           c(jb+j) = a(ia+i) - sin45*(a(ib+i)-a(id+i))
    798           d(ja+j) = -a(ic+i) - sin45*(a(ib+i)+a(id+i))
    799           d(jb+j) = a(ic+i) - sin45*(a(ib+i)+a(id+i))
    800           i = i + inc3
    801           j = j + inc4
    802        END DO
    803        ibase = ibase + inc1
    804        jbase = jbase + inc2
    805     END DO
    806 
    807     GO TO 170
    808 90  CONTINUE
    809     z = 1.0_wp/REAL(n,KIND=wp)
    810     DO l = 1, la
    811        i = ibase
    812        j = jbase
    813 !DIR$ IVDEP
    814 !CDIR NODEP
    815 !OCL NOVREC
    816        DO ijk = 1, lot
    817           c(ja+j) = z*((a(ia+i)+a(ic+i))+(a(ib+i)+a(id+i)))
    818           c(jc+j) = z*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i)))
    819           c(jb+j) = z*(a(ia+i)-a(ic+i))
    820           d(jb+j) = z*(a(id+i)-a(ib+i))
    821           i = i + inc3
    822           j = j + inc4
    823        END DO
    824        ibase = ibase + inc1
    825        jbase = jbase + inc2
    826     END DO
    827     GO TO 170
    828 
    829     ! Coding for factor 5
    830 
    831 100 CONTINUE
    832     ia = 1
    833     ib = ia + iink
    834     ic = ib + iink
    835     id = ic + iink
    836     ie = id + iink
    837     ja = 1
    838     jb = ja + (2*m-la)*inc2
    839     jc = jb + 2*m*inc2
    840     jd = jc
    841     je = jb
    842 
    843     IF (la==m) GO TO 120
    844 
    845     DO l = 1, la
    846        i = ibase
    847        j = jbase
    848 !DIR$ IVDEP
    849 !CDIR NODEP
    850 !OCL NOVREC
    851        DO ijk = 1, lot
    852           a1 = a(ib+i) + a(ie+i)
    853           a3 = a(ib+i) - a(ie+i)
    854           a2 = a(ic+i) + a(id+i)
    855           a4 = a(ic+i) - a(id+i)
    856           a5 = a(ia+i) - 0.25_wp*(a1+a2)
    857           a6 = qrt5*(a1-a2)
    858           c(ja+j) = a(ia+i) + (a1+a2)
    859           c(jb+j) = a5 + a6
    860           c(jc+j) = a5 - a6
    861           d(jb+j) = -sin72*a3 - sin36*a4
    862           d(jc+j) = -sin36*a3 + sin72*a4
    863           i = i + inc3
    864           j = j + inc4
    865        END DO
    866        ibase = ibase + inc1
    867        jbase = jbase + inc2
    868     END DO
    869     ja = ja + jink
    870     jink = 2*jink
    871     jb = jb + jink
    872     jc = jc + jink
    873     jd = jd - jink
    874     je = je - jink
    875     ibase = ibase + ijump
    876     ijump = 2*ijump + iink
    877     IF (jb==jd) GO TO 110
    878     DO k = la, kstop, la
    879        kb = k + k
    880        kc = kb + kb
    881        kd = kc + kb
    882        ke = kd + kb
    883        c1 = trigs(kb+1)
    884        s1 = trigs(kb+2)
    885        c2 = trigs(kc+1)
    886        s2 = trigs(kc+2)
    887        c3 = trigs(kd+1)
    888        s3 = trigs(kd+2)
    889        c4 = trigs(ke+1)
    890        s4 = trigs(ke+2)
    891        jbase = 0
    892        DO l = 1, la
    893           i = ibase
    894           j = jbase
    895 !DIR$ IVDEP
    896 !CDIR NODEP
    897 !OCL NOVREC
    898           DO ijk = 1, lot
    899              a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c4*a(ie+i)+s4*b(ie+i))
    900              a3 = (c1*a(ib+i)+s1*b(ib+i)) - (c4*a(ie+i)+s4*b(ie+i))
    901              a2 = (c2*a(ic+i)+s2*b(ic+i)) + (c3*a(id+i)+s3*b(id+i))
    902              a4 = (c2*a(ic+i)+s2*b(ic+i)) - (c3*a(id+i)+s3*b(id+i))
    903              b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c4*b(ie+i)-s4*a(ie+i))
    904              b3 = (c1*b(ib+i)-s1*a(ib+i)) - (c4*b(ie+i)-s4*a(ie+i))
    905              b2 = (c2*b(ic+i)-s2*a(ic+i)) + (c3*b(id+i)-s3*a(id+i))
    906              b4 = (c2*b(ic+i)-s2*a(ic+i)) - (c3*b(id+i)-s3*a(id+i))
    907              a5 = a(ia+i) - 0.25_wp*(a1+a2)
    908              a6 = qrt5*(a1-a2)
    909              b5 = b(ia+i) - 0.25_wp*(b1+b2)
    910              b6 = qrt5*(b1-b2)
    911              a10 = a5 + a6
    912              a20 = a5 - a6
    913              b10 = b5 + b6
    914              b20 = b5 - b6
    915              a11 = sin72*b3 + sin36*b4
    916              a21 = sin36*b3 - sin72*b4
    917              b11 = sin72*a3 + sin36*a4
    918              b21 = sin36*a3 - sin72*a4
    919              c(ja+j) = a(ia+i) + (a1+a2)
    920              c(jb+j) = a10 + a11
    921              c(je+j) = a10 - a11
    922              c(jc+j) = a20 + a21
    923              c(jd+j) = a20 - a21
    924              d(ja+j) = b(ia+i) + (b1+b2)
    925              d(jb+j) = b10 - b11
    926              d(je+j) = -(b10+b11)
    927              d(jc+j) = b20 - b21
    928              d(jd+j) = -(b20+b21)
    929              i = i + inc3
    930              j = j + inc4
    931           END DO
    932           ibase = ibase + inc1
    933           jbase = jbase + inc2
    934        END DO
    935        ibase = ibase + ijump
    936        ja = ja + jink
    937        jb = jb + jink
    938        jc = jc + jink
    939        jd = jd - jink
    940        je = je - jink
    941     END DO
    942     IF (jb>jd) GO TO 170
    943 110 CONTINUE
    944     jbase = 0
    945     DO l = 1, la
    946        i = ibase
    947        j = jbase
    948 !DIR$ IVDEP
    949 !CDIR NODEP
    950 !OCL NOVREC
    951        DO ijk = 1, lot
    952           a1 = a(ib+i) + a(ie+i)
    953           a3 = a(ib+i) - a(ie+i)
    954           a2 = a(ic+i) + a(id+i)
    955           a4 = a(ic+i) - a(id+i)
    956           a5 = a(ia+i) + 0.25_wp*(a3-a4)
    957           a6 = qrt5*(a3+a4)
    958           c(ja+j) = a5 + a6
    959           c(jb+j) = a5 - a6
    960           c(jc+j) = a(ia+i) - (a3-a4)
    961           d(ja+j) = -sin36*a1 - sin72*a2
    962           d(jb+j) = -sin72*a1 + sin36*a2
    963           i = i + inc3
    964           j = j + inc4
    965        END DO
    966        ibase = ibase + inc1
    967        jbase = jbase + inc2
    968     END DO
    969 
    970     GO TO 170
    971 120 CONTINUE
    972     z = 1.0_wp/REAL(n,KIND=wp)
    973     zqrt5 = z*qrt5
    974     zsin36 = z*sin36
    975     zsin72 = z*sin72
    976     DO l = 1, la
    977        i = ibase
    978        j = jbase
    979 !DIR$ IVDEP
    980 !CDIR NODEP
    981 !OCL NOVREC
    982        DO ijk = 1, lot
    983           a1 = a(ib+i) + a(ie+i)
    984           a3 = a(ib+i) - a(ie+i)
    985           a2 = a(ic+i) + a(id+i)
    986           a4 = a(ic+i) - a(id+i)
    987           a5 = z*(a(ia+i)-0.25_wp*(a1+a2))
    988           a6 = zqrt5*(a1-a2)
    989           c(ja+j) = z*(a(ia+i)+(a1+a2))
    990           c(jb+j) = a5 + a6
    991           c(jc+j) = a5 - a6
    992           d(jb+j) = -zsin72*a3 - zsin36*a4
    993           d(jc+j) = -zsin36*a3 + zsin72*a4
    994           i = i + inc3
    995           j = j + inc4
    996        END DO
    997        ibase = ibase + inc1
    998        jbase = jbase + inc2
    999     END DO
    1000     GO TO 170
    1001 
    1002     ! Coding for factor 6
    1003 
    1004 130 CONTINUE
    1005     ia = 1
    1006     ib = ia + iink
    1007     ic = ib + iink
    1008     id = ic + iink
    1009     ie = id + iink
    1010     if = ie + iink
    1011     ja = 1
    1012     jb = ja + (2*m-la)*inc2
    1013     jc = jb + 2*m*inc2
    1014     jd = jc + 2*m*inc2
    1015     je = jc
    1016     jf = jb
    1017 
    1018     IF (la==m) GO TO 150
    1019 
    1020     DO l = 1, la
    1021        i = ibase
    1022        j = jbase
    1023 !DIR$ IVDEP
    1024 !CDIR NODEP
    1025 !OCL NOVREC
    1026        DO ijk = 1, lot
    1027           a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
    1028           c(ja+j) = (a(ia+i)+a(id+i)) + a11
    1029           c(jc+j) = (a(ia+i)+a(id+i)-0.5_wp*a11)
    1030           d(jc+j) = sin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
    1031           a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
    1032           c(jb+j) = (a(ia+i)-a(id+i)) - 0.5_wp*a11
    1033           d(jb+j) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
    1034           c(jd+j) = (a(ia+i)-a(id+i)) + a11
    1035           i = i + inc3
    1036           j = j + inc4
    1037        END DO
    1038        ibase = ibase + inc1
    1039        jbase = jbase + inc2
    1040     END DO
    1041     ja = ja + jink
    1042     jink = 2*jink
    1043     jb = jb + jink
    1044     jc = jc + jink
    1045     jd = jd - jink
    1046     je = je - jink
    1047     jf = jf - jink
    1048     ibase = ibase + ijump
    1049     ijump = 2*ijump + iink
    1050     IF (jc==jd) GO TO 140
    1051     DO k = la, kstop, la
    1052        kb = k + k
    1053        kc = kb + kb
    1054        kd = kc + kb
    1055        ke = kd + kb
    1056        kf = ke + kb
    1057        c1 = trigs(kb+1)
    1058        s1 = trigs(kb+2)
    1059        c2 = trigs(kc+1)
    1060        s2 = trigs(kc+2)
    1061        c3 = trigs(kd+1)
    1062        s3 = trigs(kd+2)
    1063        c4 = trigs(ke+1)
    1064        s4 = trigs(ke+2)
    1065        c5 = trigs(kf+1)
    1066        s5 = trigs(kf+2)
    1067        jbase = 0
    1068        DO l = 1, la
    1069           i = ibase
    1070           j = jbase
    1071 !DIR$ IVDEP
    1072 !CDIR NODEP
    1073 !OCL NOVREC
    1074           DO ijk = 1, lot
    1075              a1 = c1*a(ib+i) + s1*b(ib+i)
    1076              b1 = c1*b(ib+i) - s1*a(ib+i)
    1077              a2 = c2*a(ic+i) + s2*b(ic+i)
    1078              b2 = c2*b(ic+i) - s2*a(ic+i)
    1079              a3 = c3*a(id+i) + s3*b(id+i)
    1080              b3 = c3*b(id+i) - s3*a(id+i)
    1081              a4 = c4*a(ie+i) + s4*b(ie+i)
    1082              b4 = c4*b(ie+i) - s4*a(ie+i)
    1083              a5 = c5*a(if+i) + s5*b(if+i)
    1084              b5 = c5*b(if+i) - s5*a(if+i)
    1085              a11 = (a2+a5) + (a1+a4)
    1086              a20 = (a(ia+i)+a3) - 0.5_wp*a11
    1087              a21 = sin60*((a2+a5)-(a1+a4))
    1088              b11 = (b2+b5) + (b1+b4)
    1089              b20 = (b(ia+i)+b3) - 0.5_wp*b11
    1090              b21 = sin60*((b2+b5)-(b1+b4))
    1091              c(ja+j) = (a(ia+i)+a3) + a11
    1092              d(ja+j) = (b(ia+i)+b3) + b11
    1093              c(jc+j) = a20 - b21
    1094              d(jc+j) = a21 + b20
    1095              c(je+j) = a20 + b21
    1096              d(je+j) = a21 - b20
    1097              a11 = (a2-a5) + (a4-a1)
    1098              a20 = (a(ia+i)-a3) - 0.5_wp*a11
    1099              a21 = sin60*((a4-a1)-(a2-a5))
    1100              b11 = (b5-b2) - (b4-b1)
    1101              b20 = (b3-b(ia+i)) - 0.5_wp*b11
    1102              b21 = sin60*((b5-b2)+(b4-b1))
    1103              c(jb+j) = a20 - b21
    1104              d(jb+j) = a21 - b20
    1105              c(jd+j) = a11 + (a(ia+i)-a3)
    1106              d(jd+j) = b11 + (b3-b(ia+i))
    1107              c(jf+j) = a20 + b21
    1108              d(jf+j) = a21 + b20
    1109              i = i + inc3
    1110              j = j + inc4
    1111           END DO
    1112           ibase = ibase + inc1
    1113           jbase = jbase + inc2
    1114        END DO
    1115        ibase = ibase + ijump
    1116        ja = ja + jink
    1117        jb = jb + jink
    1118        jc = jc + jink
    1119        jd = jd - jink
    1120        je = je - jink
    1121        jf = jf - jink
    1122     END DO
    1123     IF (jc>jd) GO TO 170
    1124 140 CONTINUE
    1125     jbase = 0
    1126     DO l = 1, la
    1127        i = ibase
    1128        j = jbase
    1129 !DIR$ IVDEP
    1130 !CDIR NODEP
    1131 !OCL NOVREC
    1132        DO ijk = 1, lot
    1133           c(ja+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) + sin60*(a(ib+i)-a(if+i))
    1134           d(ja+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) - sin60*(a(ic+i)+a(ie+i))
    1135           c(jb+j) = a(ia+i) - (a(ic+i)-a(ie+i))
    1136           d(jb+j) = a(id+i) - (a(ib+i)+a(if+i))
    1137           c(jc+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) - sin60*(a(ib+i)-a(if+i))
    1138           d(jc+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) + sin60*(a(ic+i)+a(ie+i))
    1139           i = i + inc3
    1140           j = j + inc4
    1141        END DO
    1142        ibase = ibase + inc1
    1143        jbase = jbase + inc2
    1144     END DO
    1145 
    1146     GO TO 170
    1147 150 CONTINUE
    1148     z = 1.0_wp/REAL(n,KIND=wp)
    1149     zsin60 = z*sin60
    1150     DO l = 1, la
    1151        i = ibase
    1152        j = jbase
    1153 !DIR$ IVDEP
    1154 !CDIR NODEP
    1155 !OCL NOVREC
    1156        DO ijk = 1, lot
    1157           a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
    1158           c(ja+j) = z*((a(ia+i)+a(id+i))+a11)
    1159           c(jc+j) = z*((a(ia+i)+a(id+i))-0.5_wp*a11)
    1160           d(jc+j) = zsin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
    1161           a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
    1162           c(jb+j) = z*((a(ia+i)-a(id+i))-0.5_wp*a11)
    1163           d(jb+j) = zsin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
    1164           c(jd+j) = z*((a(ia+i)-a(id+i))+a11)
    1165           i = i + inc3
    1166           j = j + inc4
    1167        END DO
    1168        ibase = ibase + inc1
    1169        jbase = jbase + inc2
    1170     END DO
    1171     GO TO 170
    1172 
    1173     ! Coding for factor 8
    1174 
    1175 160 CONTINUE
    1176     ibad = 3
    1177     IF (la/=m) GO TO 180
    1178     ia = 1
    1179     ib = ia + iink
    1180     ic = ib + iink
    1181     id = ic + iink
    1182     ie = id + iink
    1183     if = ie + iink
    1184     ig = if + iink
    1185     ih = ig + iink
    1186     ja = 1
    1187     jb = ja + la*inc2
    1188     jc = jb + 2*m*inc2
    1189     jd = jc + 2*m*inc2
    1190     je = jd + 2*m*inc2
    1191     z = 1.0_wp/REAL(n,KIND=wp)
    1192     zsin45 = z*SQRT(0.5_wp)
    1193 
    1194     DO l = 1, la
    1195        i = ibase
    1196        j = jbase
    1197 !DIR$ IVDEP
    1198 !CDIR NODEP
    1199 !OCL NOVREC
    1200        DO ijk = 1, lot
    1201           c(ja+j) = z*(((a(ia+i)+a(ie+i))+(a(ic+i)+a(ig+i)))+((a(id+i)+ &
    1202                &          a(ih+i))+(a(ib+i)+a(if+i))))
    1203           c(je+j) = z*(((a(ia+i)+a(ie+i))+(a(ic+i)+a(ig+i)))-((a(id+i)+ &
    1204                &          a(ih+i))+(a(ib+i)+a(if+i))))
    1205           c(jc+j) = z*((a(ia+i)+a(ie+i))-(a(ic+i)+a(ig+i)))
    1206           d(jc+j) = z*((a(id+i)+a(ih+i))-(a(ib+i)+a(if+i)))
    1207           c(jb+j) = z*(a(ia+i)-a(ie+i)) + zsin45*((a(ih+i)-a(id+i))-(a(if+ &
    1208                &          i)-a(ib+i)))
    1209           c(jd+j) = z*(a(ia+i)-a(ie+i)) - zsin45*((a(ih+i)-a(id+i))-(a(if+ &
    1210                &          i)-a(ib+i)))
    1211           d(jb+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) + &
    1212                &          z*(a(ig+i)-a(ic+i))
    1213           d(jd+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) - &
    1214                &          z*(a(ig+i)-a(ic+i))
    1215           i = i + inc3
    1216           j = j + inc4
    1217        END DO
    1218        ibase = ibase + inc1
    1219        jbase = jbase + inc2
    1220     END DO
    1221 
    1222     ! Return
    1223 
    1224 170 CONTINUE
    1225     ibad = 0
    1226 180 CONTINUE
    1227     ierr = ibad
    1228     RETURN
    1229   END SUBROUTINE qpassm
     488    IF ( igo == 7 )  igo = 6
     489
     490    IF (igo < 1 .OR. igo > 6 )  THEN
     491       ierr = 2
     492       RETURN
     493    ENDIF
     494
     495
     496    SELECT CASE ( igo )
     497
     498!
     499!--    Coding for factor 2
     500       CASE ( 1 )
     501          ia = 1
     502          ib = ia + iink
     503          ja = 1
     504          jb = ja + (2*m-la) * inc2
     505
     506          IF ( la /= m )  THEN
     507
     508             DO  l = 1, la
     509                i = ibase
     510                j = jbase
     511                !DIR$ IVDEP
     512                !CDIR NODEP
     513                !OCL NOVREC
     514                DO  ijk = 1, lot
     515                   c(ja+j) = a(ia+i) + a(ib+i)
     516                   c(jb+j) = a(ia+i) - a(ib+i)
     517                   i = i + inc3
     518                   j = j + inc4
     519                ENDDO
     520                ibase = ibase + inc1
     521                jbase = jbase + inc2
     522             ENDDO
     523
     524             ja = ja + jink
     525             jink = 2 * jink
     526             jb = jb - jink
     527             ibase = ibase + ijump
     528             ijump = 2 * ijump + iink
     529
     530             IF ( ja /= jb )  THEN
     531
     532                DO  k = la, kstop, la
     533                   kb = k + k
     534                   c1 = trigs(kb+1)
     535                   s1 = trigs(kb+2)
     536                   jbase = 0
     537                   DO  l = 1, la
     538                      i = ibase
     539                      j = jbase
     540                      !DIR$ IVDEP
     541                      !CDIR NODEP
     542                      !OCL NOVREC
     543                      DO  ijk = 1, lot
     544                         c(ja+j) = a(ia+i) + (c1*a(ib+i)+s1*b(ib+i))
     545                         c(jb+j) = a(ia+i) - (c1*a(ib+i)+s1*b(ib+i))
     546                         d(ja+j) = (c1*b(ib+i)-s1*a(ib+i)) + b(ia+i)
     547                         d(jb+j) = (c1*b(ib+i)-s1*a(ib+i)) - b(ia+i)
     548                         i = i + inc3
     549                         j = j + inc4
     550                      ENDDO
     551                      ibase = ibase + inc1
     552                      jbase = jbase + inc2
     553                   ENDDO
     554
     555                   ibase = ibase + ijump
     556                   ja = ja + jink
     557                   jb = jb - jink
     558                ENDDO
     559
     560                IF ( ja > jb )  RETURN
     561
     562             ENDIF
     563
     564             jbase = 0
     565             DO  l = 1, la
     566                i = ibase
     567                j = jbase
     568                !DIR$ IVDEP
     569                !CDIR NODEP
     570                !OCL NOVREC
     571                DO  ijk = 1, lot
     572                   c(ja+j) = a(ia+i)
     573                   d(ja+j) = -a(ib+i)
     574                   i = i + inc3
     575                   j = j + inc4
     576
     577                ENDDO
     578                ibase = ibase + inc1
     579                jbase = jbase + inc2
     580             ENDDO
     581
     582          ELSE
     583
     584             z = 1.0_wp/REAL(n,KIND=wp)
     585             DO  l = 1, la
     586                i = ibase
     587                j = jbase
     588                !DIR$ IVDEP
     589                !CDIR NODEP
     590                !OCL NOVREC
     591                DO  ijk = 1, lot
     592                   c(ja+j) = z*(a(ia+i)+a(ib+i))
     593                   c(jb+j) = z*(a(ia+i)-a(ib+i))
     594                   i = i + inc3
     595                   j = j + inc4
     596                ENDDO
     597                ibase = ibase + inc1
     598                jbase = jbase + inc2
     599             ENDDO
     600
     601          ENDIF
     602
     603!
     604!--    Coding for factor 3
     605       CASE ( 2 )
     606
     607          ia = 1
     608          ib = ia + iink
     609          ic = ib + iink
     610          ja = 1
     611          jb = ja + (2*m-la) * inc2
     612          jc = jb
     613
     614          IF ( la /= m )  THEN
     615
     616             DO  l = 1, la
     617                i = ibase
     618                j = jbase
     619                !DIR$ IVDEP
     620                !CDIR NODEP
     621                !OCL NOVREC
     622                DO  ijk = 1, lot
     623                   c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
     624                   c(jb+j) = a(ia+i) - 0.5_wp*(a(ib+i)+a(ic+i))
     625                   d(jb+j) = sin60*(a(ic+i)-a(ib+i))
     626                   i = i + inc3
     627                   j = j + inc4
     628                ENDDO
     629                ibase = ibase + inc1
     630                jbase = jbase + inc2
     631             ENDDO
     632
     633             ja = ja + jink
     634             jink = 2 * jink
     635             jb = jb + jink
     636             jc = jc - jink
     637             ibase = ibase + ijump
     638             ijump = 2 * ijump + iink
     639
     640             IF ( ja /= jc )  THEN
     641
     642                DO  k = la, kstop, la
     643                   kb = k + k
     644                   kc = kb + kb
     645                   c1 = trigs(kb+1)
     646                   s1 = trigs(kb+2)
     647                   c2 = trigs(kc+1)
     648                   s2 = trigs(kc+2)
     649
     650                   jbase = 0
     651                   DO  l = 1, la
     652                      i = ibase
     653                      j = jbase
     654                      !DIR$ IVDEP
     655                      !CDIR NODEP
     656                      !OCL NOVREC
     657                      DO  ijk = 1, lot
     658                         a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c2*a(ic+i)+s2*b(ic+i))
     659                         b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c2*b(ic+i)-s2*a(ic+i))
     660                         a2 = a(ia+i) - 0.5_wp*a1
     661                         b2 = b(ia+i) - 0.5_wp*b1
     662                         a3 = sin60*((c1*a(ib+i)+s1*b(ib+i))-(c2*a(ic+i)+s2*b(ic+i)))
     663                         b3 = sin60*((c1*b(ib+i)-s1*a(ib+i))-(c2*b(ic+i)-s2*a(ic+i)))
     664                         c(ja+j) = a(ia+i) + a1
     665                         d(ja+j) = b(ia+i) + b1
     666                         c(jb+j) = a2 + b3
     667                         d(jb+j) = b2 - a3
     668                         c(jc+j) = a2 - b3
     669                         d(jc+j) = -(b2+a3)
     670                         i = i + inc3
     671                         j = j + inc4
     672                      ENDDO
     673                      ibase = ibase + inc1
     674                      jbase = jbase + inc2
     675                   ENDDO
     676                   ibase = ibase + ijump
     677                   ja = ja + jink
     678                   jb = jb + jink
     679                   jc = jc - jink
     680                ENDDO
     681
     682                IF ( ja > jc )  RETURN
     683
     684             ENDIF
     685
     686             jbase = 0
     687             DO  l = 1, la
     688                i = ibase
     689                j = jbase
     690                !DIR$ IVDEP
     691                !CDIR NODEP
     692                !OCL NOVREC
     693                DO  ijk = 1, lot
     694                   c(ja+j) = a(ia+i) + 0.5_wp*(a(ib+i)-a(ic+i))
     695                   d(ja+j) = -sin60*(a(ib+i)+a(ic+i))
     696                   c(jb+j) = a(ia+i) - (a(ib+i)-a(ic+i))
     697                   i = i + inc3
     698                   j = j + inc4
     699                ENDDO
     700                ibase = ibase + inc1
     701                jbase = jbase + inc2
     702             ENDDO
     703
     704          ELSE
     705
     706             z = 1.0_wp / REAL( n, KIND=wp )
     707             zsin60 = z*sin60
     708             DO  l = 1, la
     709                i = ibase
     710                j = jbase
     711                !DIR$ IVDEP
     712                !CDIR NODEP
     713                !OCL NOVREC
     714                DO  ijk = 1, lot
     715                   c(ja+j) = z*(a(ia+i)+(a(ib+i)+a(ic+i)))
     716                   c(jb+j) = z*(a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))
     717                   d(jb+j) = zsin60*(a(ic+i)-a(ib+i))
     718                   i = i + inc3
     719                   j = j + inc4
     720                ENDDO
     721                ibase = ibase + inc1
     722                jbase = jbase + inc2
     723             ENDDO
     724
     725          ENDIF
     726
     727!
     728!--    Coding for factor 4
     729       CASE ( 3 )
     730
     731          ia = 1
     732          ib = ia + iink
     733          ic = ib + iink
     734          id = ic + iink
     735          ja = 1
     736          jb = ja + (2*m-la) * inc2
     737          jc = jb + 2*m*inc2
     738          jd = jb
     739
     740          IF ( la /= m )  THEN
     741
     742             DO  l = 1, la
     743
     744                i = ibase
     745                j = jbase
     746                !DIR$ IVDEP
     747                !CDIR NODEP
     748                !OCL NOVREC
     749                DO  ijk = 1, lot
     750                   c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))
     751                   c(jc+j) = (a(ia+i)+a(ic+i)) - (a(ib+i)+a(id+i))
     752                   c(jb+j) = a(ia+i) - a(ic+i)
     753                   d(jb+j) = a(id+i) - a(ib+i)
     754                   i = i + inc3
     755                   j = j + inc4
     756                ENDDO
     757                ibase = ibase + inc1
     758                jbase = jbase + inc2
     759             ENDDO
     760
     761             ja = ja + jink
     762             jink = 2 * jink
     763             jb = jb + jink
     764             jc = jc - jink
     765             jd = jd - jink
     766             ibase = ibase + ijump
     767             ijump = 2 * ijump + iink
     768
     769             IF ( jb /= jc )  THEN
     770
     771                DO  k = la, kstop, la
     772                   kb = k + k
     773                   kc = kb + kb
     774                   kd = kc + kb
     775                   c1 = trigs(kb+1)
     776                   s1 = trigs(kb+2)
     777                   c2 = trigs(kc+1)
     778                   s2 = trigs(kc+2)
     779                   c3 = trigs(kd+1)
     780                   s3 = trigs(kd+2)
     781                   jbase = 0
     782                   DO  l = 1, la
     783                      i = ibase
     784                      j = jbase
     785                      !DIR$ IVDEP
     786                      !CDIR NODEP
     787                      !OCL NOVREC
     788                      DO  ijk = 1, lot
     789                         a0 = a(ia+i) + (c2*a(ic+i)+s2*b(ic+i))
     790                         a2 = a(ia+i) - (c2*a(ic+i)+s2*b(ic+i))
     791                         a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c3*a(id+i)+s3*b(id+i))
     792                         a3 = (c1*a(ib+i)+s1*b(ib+i)) - (c3*a(id+i)+s3*b(id+i))
     793                         b0 = b(ia+i) + (c2*b(ic+i)-s2*a(ic+i))
     794                         b2 = b(ia+i) - (c2*b(ic+i)-s2*a(ic+i))
     795                         b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c3*b(id+i)-s3*a(id+i))
     796                         b3 = (c1*b(ib+i)-s1*a(ib+i)) - (c3*b(id+i)-s3*a(id+i))
     797                         c(ja+j) = a0 + a1
     798                         c(jc+j) = a0 - a1
     799                         d(ja+j) = b0 + b1
     800                         d(jc+j) = b1 - b0
     801                         c(jb+j) = a2 + b3
     802                         c(jd+j) = a2 - b3
     803                         d(jb+j) = b2 - a3
     804                         d(jd+j) = -(b2+a3)
     805                         i = i + inc3
     806                         j = j + inc4
     807                      ENDDO
     808                      ibase = ibase + inc1
     809                      jbase = jbase + inc2
     810                   ENDDO
     811                   ibase = ibase + ijump
     812                   ja = ja + jink
     813                   jb = jb + jink
     814                   jc = jc - jink
     815                   jd = jd - jink
     816                ENDDO
     817
     818                IF ( jb > jc )  RETURN
     819
     820             ENDIF
     821
     822             sin45 = SQRT( 0.5_wp )
     823             jbase = 0
     824             DO  l = 1, la
     825                i = ibase
     826                j = jbase
     827                !DIR$ IVDEP
     828                !CDIR NODEP
     829                !OCL NOVREC
     830                DO  ijk = 1, lot
     831                   c(ja+j) = a(ia+i) + sin45*(a(ib+i)-a(id+i))
     832                   c(jb+j) = a(ia+i) - sin45*(a(ib+i)-a(id+i))
     833                   d(ja+j) = -a(ic+i) - sin45*(a(ib+i)+a(id+i))
     834                   d(jb+j) = a(ic+i) - sin45*(a(ib+i)+a(id+i))
     835                   i = i + inc3
     836                   j = j + inc4
     837                ENDDO
     838                ibase = ibase + inc1
     839                jbase = jbase + inc2
     840             ENDDO
     841
     842          ELSE
     843
     844             z = 1.0_wp / REAL( n, KIND=wp )
     845             DO  l = 1, la
     846                i = ibase
     847                j = jbase
     848                !DIR$ IVDEP
     849                !CDIR NODEP
     850                !OCL NOVREC
     851                DO  ijk = 1, lot
     852                   c(ja+j) = z*((a(ia+i)+a(ic+i))+(a(ib+i)+a(id+i)))
     853                   c(jc+j) = z*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i)))
     854                   c(jb+j) = z*(a(ia+i)-a(ic+i))
     855                   d(jb+j) = z*(a(id+i)-a(ib+i))
     856                   i = i + inc3
     857                   j = j + inc4
     858                ENDDO
     859                ibase = ibase + inc1
     860                jbase = jbase + inc2
     861             ENDDO
     862
     863          ENDIF
     864
     865!
     866!--    Coding for factor 5
     867       CASE ( 4 )
     868
     869          ia = 1
     870          ib = ia + iink
     871          ic = ib + iink
     872          id = ic + iink
     873          ie = id + iink
     874          ja = 1
     875          jb = ja + (2*m-la) * inc2
     876          jc = jb + 2*m*inc2
     877          jd = jc
     878          je = jb
     879
     880          IF ( la /= m )  THEN
     881
     882             DO  l = 1, la
     883                i = ibase
     884                j = jbase
     885                !DIR$ IVDEP
     886                !CDIR NODEP
     887                !OCL NOVREC
     888                DO  ijk = 1, lot
     889                   a1 = a(ib+i) + a(ie+i)
     890                   a3 = a(ib+i) - a(ie+i)
     891                   a2 = a(ic+i) + a(id+i)
     892                   a4 = a(ic+i) - a(id+i)
     893                   a5 = a(ia+i) - 0.25_wp*(a1+a2)
     894                   a6 = qrt5*(a1-a2)
     895                   c(ja+j) = a(ia+i) + (a1+a2)
     896                   c(jb+j) = a5 + a6
     897                   c(jc+j) = a5 - a6
     898                   d(jb+j) = -sin72*a3 - sin36*a4
     899                   d(jc+j) = -sin36*a3 + sin72*a4
     900                   i = i + inc3
     901                   j = j + inc4
     902                ENDDO
     903                ibase = ibase + inc1
     904                jbase = jbase + inc2
     905             ENDDO
     906
     907             ja = ja + jink
     908             jink = 2 * jink
     909             jb = jb + jink
     910             jc = jc + jink
     911             jd = jd - jink
     912             je = je - jink
     913             ibase = ibase + ijump
     914             ijump = 2 * ijump + iink
     915
     916             IF ( jb /= jd )  THEN
     917
     918                DO  k = la, kstop, la
     919                   kb = k + k
     920                   kc = kb + kb
     921                   kd = kc + kb
     922                   ke = kd + kb
     923                   c1 = trigs(kb+1)
     924                   s1 = trigs(kb+2)
     925                   c2 = trigs(kc+1)
     926                   s2 = trigs(kc+2)
     927                   c3 = trigs(kd+1)
     928                   s3 = trigs(kd+2)
     929                   c4 = trigs(ke+1)
     930                   s4 = trigs(ke+2)
     931                   jbase = 0
     932                   DO  l = 1, la
     933                      i = ibase
     934                      j = jbase
     935                      !DIR$ IVDEP
     936                      !CDIR NODEP
     937                      !OCL NOVREC
     938                      DO  ijk = 1, lot
     939                         a1 = (c1*a(ib+i)+s1*b(ib+i)) + (c4*a(ie+i)+s4*b(ie+i))
     940                         a3 = (c1*a(ib+i)+s1*b(ib+i)) - (c4*a(ie+i)+s4*b(ie+i))
     941                         a2 = (c2*a(ic+i)+s2*b(ic+i)) + (c3*a(id+i)+s3*b(id+i))
     942                         a4 = (c2*a(ic+i)+s2*b(ic+i)) - (c3*a(id+i)+s3*b(id+i))
     943                         b1 = (c1*b(ib+i)-s1*a(ib+i)) + (c4*b(ie+i)-s4*a(ie+i))
     944                         b3 = (c1*b(ib+i)-s1*a(ib+i)) - (c4*b(ie+i)-s4*a(ie+i))
     945                         b2 = (c2*b(ic+i)-s2*a(ic+i)) + (c3*b(id+i)-s3*a(id+i))
     946                         b4 = (c2*b(ic+i)-s2*a(ic+i)) - (c3*b(id+i)-s3*a(id+i))
     947                         a5 = a(ia+i) - 0.25_wp*(a1+a2)
     948                         a6 = qrt5*(a1-a2)
     949                         b5 = b(ia+i) - 0.25_wp*(b1+b2)
     950                         b6 = qrt5*(b1-b2)
     951                         a10 = a5 + a6
     952                         a20 = a5 - a6
     953                         b10 = b5 + b6
     954                         b20 = b5 - b6
     955                         a11 = sin72*b3 + sin36*b4
     956                         a21 = sin36*b3 - sin72*b4
     957                         b11 = sin72*a3 + sin36*a4
     958                         b21 = sin36*a3 - sin72*a4
     959                         c(ja+j) = a(ia+i) + (a1+a2)
     960                         c(jb+j) = a10 + a11
     961                         c(je+j) = a10 - a11
     962                         c(jc+j) = a20 + a21
     963                         c(jd+j) = a20 - a21
     964                         d(ja+j) = b(ia+i) + (b1+b2)
     965                         d(jb+j) = b10 - b11
     966                         d(je+j) = -(b10+b11)
     967                         d(jc+j) = b20 - b21
     968                         d(jd+j) = -(b20+b21)
     969                         i = i + inc3
     970                         j = j + inc4
     971                      ENDDO
     972                      ibase = ibase + inc1
     973                      jbase = jbase + inc2
     974                   ENDDO
     975                   ibase = ibase + ijump
     976                   ja = ja + jink
     977                   jb = jb + jink
     978                   jc = jc + jink
     979                   jd = jd - jink
     980                   je = je - jink
     981                ENDDO
     982
     983                IF ( jb > jd )  RETURN
     984
     985             ENDIF
     986
     987             jbase = 0
     988             DO  l = 1, la
     989                i = ibase
     990                j = jbase
     991                !DIR$ IVDEP
     992                !CDIR NODEP
     993                !OCL NOVREC
     994                DO  ijk = 1, lot
     995                   a1 = a(ib+i) + a(ie+i)
     996                   a3 = a(ib+i) - a(ie+i)
     997                   a2 = a(ic+i) + a(id+i)
     998                   a4 = a(ic+i) - a(id+i)
     999                   a5 = a(ia+i) + 0.25_wp*(a3-a4)
     1000                   a6 = qrt5*(a3+a4)
     1001                   c(ja+j) = a5 + a6
     1002                   c(jb+j) = a5 - a6
     1003                   c(jc+j) = a(ia+i) - (a3-a4)
     1004                   d(ja+j) = -sin36*a1 - sin72*a2
     1005                   d(jb+j) = -sin72*a1 + sin36*a2
     1006                   i = i + inc3
     1007                   j = j + inc4
     1008                ENDDO
     1009                ibase = ibase + inc1
     1010                jbase = jbase + inc2
     1011             ENDDO
     1012
     1013          ELSE
     1014
     1015             z = 1.0_wp / REAL( n, KIND=wp )
     1016             zqrt5  = z * qrt5
     1017             zsin36 = z * sin36
     1018             zsin72 = z * sin72
     1019             DO  l = 1, la
     1020                i = ibase
     1021                j = jbase
     1022                !DIR$ IVDEP
     1023                !CDIR NODEP
     1024                !OCL NOVREC
     1025                DO  ijk = 1, lot
     1026                   a1 = a(ib+i) + a(ie+i)
     1027                   a3 = a(ib+i) - a(ie+i)
     1028                   a2 = a(ic+i) + a(id+i)
     1029                   a4 = a(ic+i) - a(id+i)
     1030                   a5 = z*(a(ia+i)-0.25_wp*(a1+a2))
     1031                   a6 = zqrt5*(a1-a2)
     1032                   c(ja+j) = z*(a(ia+i)+(a1+a2))
     1033                   c(jb+j) = a5 + a6
     1034                   c(jc+j) = a5 - a6
     1035                   d(jb+j) = -zsin72*a3 - zsin36*a4
     1036                   d(jc+j) = -zsin36*a3 + zsin72*a4
     1037                   i = i + inc3
     1038                   j = j + inc4
     1039                ENDDO
     1040                ibase = ibase + inc1
     1041                jbase = jbase + inc2
     1042             ENDDO
     1043
     1044          ENDIF
     1045
     1046!
     1047!--    Coding for factor 6
     1048       CASE ( 5 )
     1049
     1050          ia = 1
     1051          ib = ia + iink
     1052          ic = ib + iink
     1053          id = ic + iink
     1054          ie = id + iink
     1055          if = ie + iink
     1056          ja = 1
     1057          jb = ja + (2*m-la) * inc2
     1058          jc = jb + 2*m*inc2
     1059          jd = jc + 2*m*inc2
     1060          je = jc
     1061          jf = jb
     1062
     1063          IF ( la /= m )  THEN
     1064
     1065             DO  l = 1, la
     1066                i = ibase
     1067                j = jbase
     1068                !DIR$ IVDEP
     1069                !CDIR NODEP
     1070                !OCL NOVREC
     1071                DO  ijk = 1, lot
     1072                   a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
     1073                   c(ja+j) = (a(ia+i)+a(id+i)) + a11
     1074                   c(jc+j) = (a(ia+i)+a(id+i)-0.5_wp*a11)
     1075                   d(jc+j) = sin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
     1076                   a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
     1077                   c(jb+j) = (a(ia+i)-a(id+i)) - 0.5_wp*a11
     1078                   d(jb+j) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
     1079                   c(jd+j) = (a(ia+i)-a(id+i)) + a11
     1080                   i = i + inc3
     1081                   j = j + inc4
     1082                ENDDO
     1083                ibase = ibase + inc1
     1084                jbase = jbase + inc2
     1085             ENDDO
     1086             ja = ja + jink
     1087             jink = 2 * jink
     1088             jb = jb + jink
     1089             jc = jc + jink
     1090             jd = jd - jink
     1091             je = je - jink
     1092             jf = jf - jink
     1093             ibase = ibase + ijump
     1094             ijump = 2 * ijump + iink
     1095
     1096             IF ( jc /= jd )  THEN
     1097
     1098                DO  k = la, kstop, la
     1099                   kb = k + k
     1100                   kc = kb + kb
     1101                   kd = kc + kb
     1102                   ke = kd + kb
     1103                   kf = ke + kb
     1104                   c1 = trigs(kb+1)
     1105                   s1 = trigs(kb+2)
     1106                   c2 = trigs(kc+1)
     1107                   s2 = trigs(kc+2)
     1108                   c3 = trigs(kd+1)
     1109                   s3 = trigs(kd+2)
     1110                   c4 = trigs(ke+1)
     1111                   s4 = trigs(ke+2)
     1112                   c5 = trigs(kf+1)
     1113                   s5 = trigs(kf+2)
     1114                   jbase = 0
     1115                   DO  l = 1, la
     1116                      i = ibase
     1117                      j = jbase
     1118                      !DIR$ IVDEP
     1119                      !CDIR NODEP
     1120                      !OCL NOVREC
     1121                      DO  ijk = 1, lot
     1122                         a1 = c1*a(ib+i) + s1*b(ib+i)
     1123                         b1 = c1*b(ib+i) - s1*a(ib+i)
     1124                         a2 = c2*a(ic+i) + s2*b(ic+i)
     1125                         b2 = c2*b(ic+i) - s2*a(ic+i)
     1126                         a3 = c3*a(id+i) + s3*b(id+i)
     1127                         b3 = c3*b(id+i) - s3*a(id+i)
     1128                         a4 = c4*a(ie+i) + s4*b(ie+i)
     1129                         b4 = c4*b(ie+i) - s4*a(ie+i)
     1130                         a5 = c5*a(if+i) + s5*b(if+i)
     1131                         b5 = c5*b(if+i) - s5*a(if+i)
     1132                         a11 = (a2+a5) + (a1+a4)
     1133                         a20 = (a(ia+i)+a3) - 0.5_wp*a11
     1134                         a21 = sin60*((a2+a5)-(a1+a4))
     1135                         b11 = (b2+b5) + (b1+b4)
     1136                         b20 = (b(ia+i)+b3) - 0.5_wp*b11
     1137                         b21 = sin60*((b2+b5)-(b1+b4))
     1138                         c(ja+j) = (a(ia+i)+a3) + a11
     1139                         d(ja+j) = (b(ia+i)+b3) + b11
     1140                         c(jc+j) = a20 - b21
     1141                         d(jc+j) = a21 + b20
     1142                         c(je+j) = a20 + b21
     1143                         d(je+j) = a21 - b20
     1144                         a11 = (a2-a5) + (a4-a1)
     1145                         a20 = (a(ia+i)-a3) - 0.5_wp*a11
     1146                         a21 = sin60*((a4-a1)-(a2-a5))
     1147                         b11 = (b5-b2) - (b4-b1)
     1148                         b20 = (b3-b(ia+i)) - 0.5_wp*b11
     1149                         b21 = sin60*((b5-b2)+(b4-b1))
     1150                         c(jb+j) = a20 - b21
     1151                         d(jb+j) = a21 - b20
     1152                         c(jd+j) = a11 + (a(ia+i)-a3)
     1153                         d(jd+j) = b11 + (b3-b(ia+i))
     1154                         c(jf+j) = a20 + b21
     1155                         d(jf+j) = a21 + b20
     1156                         i = i + inc3
     1157                         j = j + inc4
     1158                      ENDDO
     1159                      ibase = ibase + inc1
     1160                      jbase = jbase + inc2
     1161                   ENDDO
     1162                   ibase = ibase + ijump
     1163                   ja = ja + jink
     1164                   jb = jb + jink
     1165                   jc = jc + jink
     1166                   jd = jd - jink
     1167                   je = je - jink
     1168                   jf = jf - jink
     1169                ENDDO
     1170
     1171                IF ( jc > jd )  RETURN
     1172
     1173             ENDIF
     1174
     1175             jbase = 0
     1176             DO  l = 1, la
     1177                i = ibase
     1178                j = jbase
     1179                !DIR$ IVDEP
     1180                !CDIR NODEP
     1181                !OCL NOVREC
     1182                DO  ijk = 1, lot
     1183                   c(ja+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) + sin60*(a(ib+i)-a(if+i))
     1184                   d(ja+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) - sin60*(a(ic+i)+a(ie+i))
     1185                   c(jb+j) = a(ia+i) - (a(ic+i)-a(ie+i))
     1186                   d(jb+j) = a(id+i) - (a(ib+i)+a(if+i))
     1187                   c(jc+j) = (a(ia+i)+0.5_wp*(a(ic+i)-a(ie+i))) - sin60*(a(ib+i)-a(if+i))
     1188                   d(jc+j) = -(a(id+i)+0.5_wp*(a(ib+i)+a(if+i))) + sin60*(a(ic+i)+a(ie+i))
     1189                   i = i + inc3
     1190                   j = j + inc4
     1191                ENDDO
     1192                ibase = ibase + inc1
     1193                jbase = jbase + inc2
     1194             ENDDO
     1195
     1196          ELSE
     1197
     1198             z = 1.0_wp/REAL(n,KIND=wp)
     1199             zsin60 = z*sin60
     1200             DO  l = 1, la
     1201                i = ibase
     1202                j = jbase
     1203                !DIR$ IVDEP
     1204                !CDIR NODEP
     1205                !OCL NOVREC
     1206                DO  ijk = 1, lot
     1207                   a11 = (a(ic+i)+a(if+i)) + (a(ib+i)+a(ie+i))
     1208                   c(ja+j) = z*((a(ia+i)+a(id+i))+a11)
     1209                   c(jc+j) = z*((a(ia+i)+a(id+i))-0.5_wp*a11)
     1210                   d(jc+j) = zsin60*((a(ic+i)+a(if+i))-(a(ib+i)+a(ie+i)))
     1211                   a11 = (a(ic+i)-a(if+i)) + (a(ie+i)-a(ib+i))
     1212                   c(jb+j) = z*((a(ia+i)-a(id+i))-0.5_wp*a11)
     1213                   d(jb+j) = zsin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
     1214                   c(jd+j) = z*((a(ia+i)-a(id+i))+a11)
     1215                   i = i + inc3
     1216                   j = j + inc4
     1217                ENDDO
     1218                ibase = ibase + inc1
     1219                jbase = jbase + inc2
     1220             ENDDO
     1221
     1222          ENDIF
     1223
     1224!
     1225!--    Coding for factor 8
     1226       CASE ( 6 )
     1227
     1228          IF ( la /= m )  THEN
     1229             ierr = 3
     1230             RETURN
     1231          ENDIF
     1232
     1233          ia = 1
     1234          ib = ia + iink
     1235          ic = ib + iink
     1236          id = ic + iink
     1237          ie = id + iink
     1238          if = ie + iink
     1239          ig = if + iink
     1240          ih = ig + iink
     1241          ja = 1
     1242          jb = ja + la * inc2
     1243          jc = jb + 2*m*inc2
     1244          jd = jc + 2*m*inc2
     1245          je = jd + 2*m*inc2
     1246          z = 1.0_wp / REAL( n, KIND=wp )
     1247          zsin45 = z * SQRT( 0.5_wp )
     1248
     1249          DO  l = 1, la
     1250             i = ibase
     1251             j = jbase
     1252
     1253             !DIR$ IVDEP
     1254             !CDIR NODEP
     1255             !OCL NOVREC
     1256             DO  ijk = 1, lot
     1257                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))))
     1258                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))))
     1259                c(jc+j) = z*((a(ia+i)+a(ie+i))-(a(ic+i)+a(ig+i)))
     1260                d(jc+j) = z*((a(id+i)+a(ih+i))-(a(ib+i)+a(if+i)))
     1261                c(jb+j) = z*(a(ia+i)-a(ie+i)) + zsin45*((a(ih+i)-a(id+i))-(a(if+i)-a(ib+i)))
     1262                c(jd+j) = z*(a(ia+i)-a(ie+i)) - zsin45*((a(ih+i)-a(id+i))-(a(if+i)-a(ib+i)))
     1263                d(jb+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) + z*(a(ig+i)-a(ic+i))
     1264                d(jd+j) = zsin45*((a(ih+i)-a(id+i))+(a(if+i)-a(ib+i))) - z*(a(ig+i)-a(ic+i))
     1265                i = i + inc3
     1266                j = j + inc4
     1267             ENDDO
     1268             ibase = ibase + inc1
     1269             jbase = jbase + inc2
     1270          ENDDO
     1271
     1272    END SELECT
     1273
     1274 END SUBROUTINE qpassm
    12301275
    12311276!------------------------------------------------------------------------------!
    12321277! Description:
    12331278! ------------
    1234 !> @todo Missing subroutine description.
     1279!> Same as qpassm, but for backward fft
    12351280!------------------------------------------------------------------------------!
    1236   SUBROUTINE rpassm(a,b,c,d,trigs,inc1,inc2,inc3,inc4,lot,n,ifac,la,ierr)
    1237     ! Dimension a(n),b(n),c(n),d(n),trigs(n)
     1281 SUBROUTINE rpassm( a, b, c, d, trigs, inc1, inc2, inc3, inc4, lot, n, ifac, la, ierr )
    12381282
    12391283    USE kinds
     
    12411285    IMPLICIT NONE
    12421286
    1243     !  Scalar arguments
    12441287    INTEGER(iwp) ::  ierr !<
    12451288    INTEGER(iwp) ::  ifac !<
     
    12511294    INTEGER(iwp) ::  lot  !<
    12521295    INTEGER(iwp) ::  n    !<
    1253 
    1254     !  Array arguments
     1296!
     1297!-- Arrays are dimensioned with n
    12551298    REAL(wp) ::  a(*)     !<
    12561299    REAL(wp) ::  b(*)     !<
     
    12591302    REAL(wp) ::  trigs(*) !<
    12601303
    1261     !  Local scalars:
    12621304    REAL(wp) ::  c1     !<
    12631305    REAL(wp) ::  c2     !<
     
    12841326    INTEGER(iwp) ::  ia    !<
    12851327    INTEGER(iwp) ::  ib    !<
    1286     INTEGER(iwp) ::  ibad  !<
    12871328    INTEGER(iwp) ::  ibase !<
    12881329    INTEGER(iwp) ::  ic    !<
     
    13151356    INTEGER(iwp) ::  m     !<
    13161357
    1317     !  Local arrays:
    13181358    REAL(wp) ::  a10(nfft) !<
    13191359    REAL(wp) ::  a11(nfft) !<
     
    13251365    REAL(wp) ::  b21(nfft) !<
    13261366
    1327     !  Intrinsic functions
    1328 !    INTRINSIC SQRT
    1329 
    1330     !  Data statements
    1331     DATA sin36/0.587785252292473_wp/, sin72/0.951056516295154_wp/, &
    1332          &      qrt5/0.559016994374947_wp/, sin60/0.866025403784437_wp/
    1333 
    1334 
    1335     !  Executable statements
    1336 
    1337     m = n/ifac
    1338     iink = la*inc1
    1339     jink = la*inc2
    1340     jump = (ifac-1)*jink
    1341     kstop = (n-ifac)/(2*ifac)
    1342 
    1343     ibad = 1
    1344     IF (lot>nfft) GO TO 180
     1367
     1368    DATA  sin36/0.587785252292473_wp/, sin72/0.951056516295154_wp/,                                &
     1369          qrt5/0.559016994374947_wp/,  sin60/0.866025403784437_wp/
     1370
     1371
     1372    ierr = 0
     1373
     1374    m = n / ifac
     1375    iink = la * inc1
     1376    jink = la * inc2
     1377    jump = (ifac-1) * jink
     1378    kstop = (n-ifac) / (2*ifac)
     1379
     1380    IF ( lot > nfft )  THEN
     1381       ierr = 1
     1382       RETURN
     1383    ENDIF
    13451384    ibase = 0
    13461385    jbase = 0
    13471386    igo = ifac - 1
    1348     IF (igo==7) igo = 6
    1349     ibad = 2
    1350     IF (igo<1 .OR. igo>6) GO TO 180
    1351     GO TO (10,40,70,100,130,160) igo
    1352 
    1353     ! Coding for factor 2
    1354 
    1355 10  CONTINUE
    1356     ia = 1
    1357     ib = ia + (2*m-la)*inc1
    1358     ja = 1
    1359     jb = ja + jink
    1360 
    1361     IF (la==m) GO TO 30
    1362 
    1363     DO l = 1, la
    1364        i = ibase
    1365        j = jbase
    1366 !DIR$ IVDEP
    1367 !CDIR NODEP
    1368 !OCL NOVREC
    1369        DO ijk = 1, lot
    1370           c(ja+j) = a(ia+i) + a(ib+i)
    1371           c(jb+j) = a(ia+i) - a(ib+i)
    1372           i = i + inc3
    1373           j = j + inc4
    1374        END DO
    1375        ibase = ibase + inc1
    1376        jbase = jbase + inc2
    1377     END DO
    1378     ia = ia + iink
    1379     iink = 2*iink
    1380     ib = ib - iink
    1381     ibase = 0
    1382     jbase = jbase + jump
    1383     jump = 2*jump + jink
    1384     IF (ia==ib) GO TO 20
    1385     DO k = la, kstop, la
    1386        kb = k + k
    1387        c1 = trigs(kb+1)
    1388        s1 = trigs(kb+2)
    1389        ibase = 0
    1390        DO l = 1, la
    1391           i = ibase
    1392           j = jbase
    1393 !DIR$ IVDEP
    1394 !CDIR NODEP
    1395 !OCL NOVREC
    1396           DO ijk = 1, lot
    1397              c(ja+j) = a(ia+i) + a(ib+i)
    1398              d(ja+j) = b(ia+i) - b(ib+i)
    1399              c(jb+j) = c1*(a(ia+i)-a(ib+i)) - s1*(b(ia+i)+b(ib+i))
    1400              d(jb+j) = s1*(a(ia+i)-a(ib+i)) + c1*(b(ia+i)+b(ib+i))
    1401              i = i + inc3
    1402              j = j + inc4
    1403           END DO
    1404           ibase = ibase + inc1
    1405           jbase = jbase + inc2
    1406        END DO
    1407        ia = ia + iink
    1408        ib = ib - iink
    1409        jbase = jbase + jump
    1410     END DO
    1411     IF (ia>ib) GO TO 170
    1412 20  CONTINUE
    1413     ibase = 0
    1414     DO l = 1, la
    1415        i = ibase
    1416        j = jbase
    1417 !DIR$ IVDEP
    1418 !CDIR NODEP
    1419 !OCL NOVREC
    1420        DO ijk = 1, lot
    1421           c(ja+j) = a(ia+i)
    1422           c(jb+j) = -b(ia+i)
    1423           i = i + inc3
    1424           j = j + inc4
    1425        END DO
    1426        ibase = ibase + inc1
    1427        jbase = jbase + inc2
    1428     END DO
    1429 
    1430     GO TO 170
    1431 30  CONTINUE
    1432     DO l = 1, la
    1433        i = ibase
    1434        j = jbase
    1435 !DIR$ IVDEP
    1436 !CDIR NODEP
    1437 !OCL NOVREC
    1438        DO ijk = 1, lot
    1439           c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i))
    1440           c(jb+j) = 2.0_wp*(a(ia+i)-a(ib+i))
    1441           i = i + inc3
    1442           j = j + inc4
    1443        END DO
    1444        ibase = ibase + inc1
    1445        jbase = jbase + inc2
    1446     END DO
    1447     GO TO 170
    1448 
    1449     ! Coding for factor 3
    1450 
    1451 40  CONTINUE
    1452     ia = 1
    1453     ib = ia + (2*m-la)*inc1
    1454     ic = ib
    1455     ja = 1
    1456     jb = ja + jink
    1457     jc = jb + jink
    1458 
    1459     IF (la==m) GO TO 60
    1460 
    1461     DO l = 1, la
    1462        i = ibase
    1463        j = jbase
    1464 !DIR$ IVDEP
    1465 !CDIR NODEP
    1466 !OCL NOVREC
    1467        DO ijk = 1, lot
    1468           c(ja+j) = a(ia+i) + a(ib+i)
    1469           c(jb+j) = (a(ia+i)-0.5_wp*a(ib+i)) - (sin60*(b(ib+i)))
    1470           c(jc+j) = (a(ia+i)-0.5_wp*a(ib+i)) + (sin60*(b(ib+i)))
    1471           i = i + inc3
    1472           j = j + inc4
    1473        END DO
    1474        ibase = ibase + inc1
    1475        jbase = jbase + inc2
    1476     END DO
    1477     ia = ia + iink
    1478     iink = 2*iink
    1479     ib = ib + iink
    1480     ic = ic - iink
    1481     jbase = jbase + jump
    1482     jump = 2*jump + jink
    1483     IF (ia==ic) GO TO 50
    1484     DO k = la, kstop, la
    1485        kb = k + k
    1486        kc = kb + kb
    1487        c1 = trigs(kb+1)
    1488        s1 = trigs(kb+2)
    1489        c2 = trigs(kc+1)
    1490        s2 = trigs(kc+2)
    1491        ibase = 0
    1492        DO l = 1, la
    1493           i = ibase
    1494           j = jbase
    1495 !DIR$ IVDEP
    1496 !CDIR NODEP
    1497 !OCL NOVREC
    1498           DO ijk = 1, lot
    1499              c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
    1500              d(ja+j) = b(ia+i) + (b(ib+i)-b(ic+i))
    1501              c(jb+j) = c1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ &
    1502                   &            b(ic+i))))                                      &
    1503                   &    - s1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- &
    1504                   &            a(ic+i))))
    1505              d(jb+j) = s1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ &
    1506                   &            b(ic+i))))                                      &
    1507                   &    + c1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- &
    1508                   &            a(ic+i))))
    1509              c(jc+j) = c2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ &
    1510                   &            b(ic+i))))                                      &
    1511                   &    - s2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- &
    1512                   &            a(ic+i))))
    1513              d(jc+j) = s2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ &
    1514                   &            b(ic+i))))                                      &
    1515                   &    + c2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- &
    1516                   &            a(ic+i))))
    1517              i = i + inc3
    1518              j = j + inc4
    1519           END DO
    1520           ibase = ibase + inc1
    1521           jbase = jbase + inc2
    1522        END DO
    1523        ia = ia + iink
    1524        ib = ib + iink
    1525        ic = ic - iink
    1526        jbase = jbase + jump
    1527     END DO
    1528     IF (ia>ic) GO TO 170
    1529 50  CONTINUE
    1530     ibase = 0
    1531     DO l = 1, la
    1532        i = ibase
    1533        j = jbase
    1534 !DIR$ IVDEP
    1535 !CDIR NODEP
    1536 !OCL NOVREC
    1537        DO ijk = 1, lot
    1538           c(ja+j) = a(ia+i) + a(ib+i)
    1539           c(jb+j) = (0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
    1540           c(jc+j) = -(0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
    1541           i = i + inc3
    1542           j = j + inc4
    1543        END DO
    1544        ibase = ibase + inc1
    1545        jbase = jbase + inc2
    1546     END DO
    1547 
    1548     GO TO 170
    1549 60  CONTINUE
    1550     ssin60 = 2.0_wp*sin60
    1551     DO l = 1, la
    1552        i = ibase
    1553        j = jbase
    1554 !DIR$ IVDEP
    1555 !CDIR NODEP
    1556 !OCL NOVREC
    1557        DO ijk = 1, lot
    1558           c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i))
    1559           c(jb+j) = (2.0_wp*a(ia+i)-a(ib+i)) - (ssin60*b(ib+i))
    1560           c(jc+j) = (2.0_wp*a(ia+i)-a(ib+i)) + (ssin60*b(ib+i))
    1561           i = i + inc3
    1562           j = j + inc4
    1563        END DO
    1564        ibase = ibase + inc1
    1565        jbase = jbase + inc2
    1566     END DO
    1567     GO TO 170
    1568 
    1569     ! Coding for factor 4
    1570 
    1571 70  CONTINUE
    1572     ia = 1
    1573     ib = ia + (2*m-la)*inc1
    1574     ic = ib + 2*m*inc1
    1575     id = ib
    1576     ja = 1
    1577     jb = ja + jink
    1578     jc = jb + jink
    1579     jd = jc + jink
    1580 
    1581     IF (la==m) GO TO 90
    1582 
    1583     DO l = 1, la
    1584        i = ibase
    1585        j = jbase
    1586 !DIR$ IVDEP
    1587 !CDIR NODEP
    1588 !OCL NOVREC
    1589        DO ijk = 1, lot
    1590           c(ja+j) = (a(ia+i)+a(ic+i)) + a(ib+i)
    1591           c(jb+j) = (a(ia+i)-a(ic+i)) - b(ib+i)
    1592           c(jc+j) = (a(ia+i)+a(ic+i)) - a(ib+i)
    1593           c(jd+j) = (a(ia+i)-a(ic+i)) + b(ib+i)
    1594           i = i + inc3
    1595           j = j + inc4
    1596        END DO
    1597        ibase = ibase + inc1
    1598        jbase = jbase + inc2
    1599     END DO
    1600     ia = ia + iink
    1601     iink = 2*iink
    1602     ib = ib + iink
    1603     ic = ic - iink
    1604     id = id - iink
    1605     jbase = jbase + jump
    1606     jump = 2*jump + jink
    1607     IF (ib==ic) GO TO 80
    1608     DO k = la, kstop, la
    1609        kb = k + k
    1610        kc = kb + kb
    1611        kd = kc + kb
    1612        c1 = trigs(kb+1)
    1613        s1 = trigs(kb+2)
    1614        c2 = trigs(kc+1)
    1615        s2 = trigs(kc+2)
    1616        c3 = trigs(kd+1)
    1617        s3 = trigs(kd+2)
    1618        ibase = 0
    1619        DO l = 1, la
    1620           i = ibase
    1621           j = jbase
    1622 !DIR$ IVDEP
    1623 !CDIR NODEP
    1624 !OCL NOVREC
    1625           DO ijk = 1, lot
    1626              c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))
    1627              d(ja+j) = (b(ia+i)-b(ic+i)) + (b(ib+i)-b(id+i))
    1628              c(jc+j) = c2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) - s2*((b(ia+ &
    1629                   &            i)-b(ic+i))-(b(ib+i)-b(id+i)))
    1630              d(jc+j) = s2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) + c2*((b(ia+ &
    1631                   &            i)-b(ic+i))-(b(ib+i)-b(id+i)))
    1632              c(jb+j) = c1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) - s1*((b(ia+ &
    1633                   &            i)+b(ic+i))+(a(ib+i)-a(id+i)))
    1634              d(jb+j) = s1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) + c1*((b(ia+ &
    1635                   &            i)+b(ic+i))+(a(ib+i)-a(id+i)))
    1636              c(jd+j) = c3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) - s3*((b(ia+ &
    1637                   &            i)+b(ic+i))-(a(ib+i)-a(id+i)))
    1638              d(jd+j) = s3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) + c3*((b(ia+ &
    1639                   &            i)+b(ic+i))-(a(ib+i)-a(id+i)))
    1640              i = i + inc3
    1641              j = j + inc4
    1642           END DO
    1643           ibase = ibase + inc1
    1644           jbase = jbase + inc2
    1645        END DO
    1646        ia = ia + iink
    1647        ib = ib + iink
    1648        ic = ic - iink
    1649        id = id - iink
    1650        jbase = jbase + jump
    1651     END DO
    1652     IF (ib>ic) GO TO 170
    1653 80  CONTINUE
    1654     ibase = 0
    1655     sin45 = SQRT(0.5_wp)
    1656     DO l = 1, la
    1657        i = ibase
    1658        j = jbase
    1659 !DIR$ IVDEP
    1660 !CDIR NODEP
    1661 !OCL NOVREC
    1662        DO ijk = 1, lot
    1663           c(ja+j) = a(ia+i) + a(ib+i)
    1664           c(jb+j) = sin45*((a(ia+i)-a(ib+i))-(b(ia+i)+b(ib+i)))
    1665           c(jc+j) = b(ib+i) - b(ia+i)
    1666           c(jd+j) = -sin45*((a(ia+i)-a(ib+i))+(b(ia+i)+b(ib+i)))
    1667           i = i + inc3
    1668           j = j + inc4
    1669        END DO
    1670        ibase = ibase + inc1
    1671        jbase = jbase + inc2
    1672     END DO
    1673 
    1674     GO TO 170
    1675 90  CONTINUE
    1676     DO l = 1, la
    1677        i = ibase
    1678        j = jbase
    1679 !DIR$ IVDEP
    1680 !CDIR NODEP
    1681 !OCL NOVREC
    1682        DO ijk = 1, lot
    1683           c(ja+j) = 2.0_wp*((a(ia+i)+a(ic+i))+a(ib+i))
    1684           c(jb+j) = 2.0_wp*((a(ia+i)-a(ic+i))-b(ib+i))
    1685           c(jc+j) = 2.0_wp*((a(ia+i)+a(ic+i))-a(ib+i))
    1686           c(jd+j) = 2.0_wp*((a(ia+i)-a(ic+i))+b(ib+i))
    1687           i = i + inc3
    1688           j = j + inc4
    1689        END DO
    1690        ibase = ibase + inc1
    1691        jbase = jbase + inc2
    1692     END DO
    1693 
    1694     ! Coding for factor 5
    1695 
    1696     GO TO 170
    1697 100 CONTINUE
    1698     ia = 1
    1699     ib = ia + (2*m-la)*inc1
    1700     ic = ib + 2*m*inc1
    1701     id = ic
    1702     ie = ib
    1703     ja = 1
    1704     jb = ja + jink
    1705     jc = jb + jink
    1706     jd = jc + jink
    1707     je = jd + jink
    1708 
    1709     IF (la==m) GO TO 120
    1710 
    1711     DO l = 1, la
    1712        i = ibase
    1713        j = jbase
    1714 !DIR$ IVDEP
    1715 !CDIR NODEP
    1716 !OCL NOVREC
    1717        DO ijk = 1, lot
    1718           c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
    1719           c(jb+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) - &
    1720                &          (sin72*b(ib+i)+sin36*b(ic+i))
    1721           c(jc+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) - &
    1722                &          (sin36*b(ib+i)-sin72*b(ic+i))
    1723           c(jd+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) + &
    1724                &          (sin36*b(ib+i)-sin72*b(ic+i))
    1725           c(je+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) + &
    1726                &          (sin72*b(ib+i)+sin36*b(ic+i))
    1727           i = i + inc3
    1728           j = j + inc4
    1729        END DO
    1730        ibase = ibase + inc1
    1731        jbase = jbase + inc2
    1732     END DO
    1733     ia = ia + iink
    1734     iink = 2*iink
    1735     ib = ib + iink
    1736     ic = ic + iink
    1737     id = id - iink
    1738     ie = ie - iink
    1739     jbase = jbase + jump
    1740     jump = 2*jump + jink
    1741     IF (ib==id) GO TO 110
    1742     DO k = la, kstop, la
    1743        kb = k + k
    1744        kc = kb + kb
    1745        kd = kc + kb
    1746        ke = kd + kb
    1747        c1 = trigs(kb+1)
    1748        s1 = trigs(kb+2)
    1749        c2 = trigs(kc+1)
    1750        s2 = trigs(kc+2)
    1751        c3 = trigs(kd+1)
    1752        s3 = trigs(kd+2)
    1753        c4 = trigs(ke+1)
    1754        s4 = trigs(ke+2)
    1755        ibase = 0
    1756        DO l = 1, la
    1757           i = ibase
    1758           j = jbase
    1759 !DIR$ IVDEP
    1760 !CDIR NODEP
    1761 !OCL NOVREC
    1762           DO ijk = 1, lot
    1763 
    1764              a10(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) + &
    1765                   &            qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
    1766              a20(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) - &
    1767                   &            qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
    1768              b10(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) + &
    1769                   &            qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
    1770              b20(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) - &
    1771                   &            qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
    1772              a11(ijk) = sin72*(b(ib+i)+b(ie+i)) + sin36*(b(ic+i)+b(id+i))
    1773              a21(ijk) = sin36*(b(ib+i)+b(ie+i)) - sin72*(b(ic+i)+b(id+i))
    1774              b11(ijk) = sin72*(a(ib+i)-a(ie+i)) + sin36*(a(ic+i)-a(id+i))
    1775              b21(ijk) = sin36*(a(ib+i)-a(ie+i)) - sin72*(a(ic+i)-a(id+i))
    1776 
    1777              c(ja+j) = a(ia+i) + ((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))
    1778              d(ja+j) = b(ia+i) + ((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))
    1779              c(jb+j) = c1*(a10(ijk)-a11(ijk)) - s1*(b10(ijk)+b11(ijk))
    1780              d(jb+j) = s1*(a10(ijk)-a11(ijk)) + c1*(b10(ijk)+b11(ijk))
    1781              c(je+j) = c4*(a10(ijk)+a11(ijk)) - s4*(b10(ijk)-b11(ijk))
    1782              d(je+j) = s4*(a10(ijk)+a11(ijk)) + c4*(b10(ijk)-b11(ijk))
    1783              c(jc+j) = c2*(a20(ijk)-a21(ijk)) - s2*(b20(ijk)+b21(ijk))
    1784              d(jc+j) = s2*(a20(ijk)-a21(ijk)) + c2*(b20(ijk)+b21(ijk))
    1785              c(jd+j) = c3*(a20(ijk)+a21(ijk)) - s3*(b20(ijk)-b21(ijk))
    1786              d(jd+j) = s3*(a20(ijk)+a21(ijk)) + c3*(b20(ijk)-b21(ijk))
    1787 
    1788              i = i + inc3
    1789              j = j + inc4
    1790           END DO
    1791           ibase = ibase + inc1
    1792           jbase = jbase + inc2
    1793        END DO
    1794        ia = ia + iink
    1795        ib = ib + iink
    1796        ic = ic + iink
    1797        id = id - iink
    1798        ie = ie - iink
    1799        jbase = jbase + jump
    1800     END DO
    1801     IF (ib>id) GO TO 170
    1802 110 CONTINUE
    1803     ibase = 0
    1804     DO l = 1, la
    1805        i = ibase
    1806        j = jbase
    1807 !DIR$ IVDEP
    1808 !CDIR NODEP
    1809 !OCL NOVREC
    1810        DO ijk = 1, lot
    1811           c(ja+j) = (a(ia+i)+a(ib+i)) + a(ic+i)
    1812           c(jb+j) = (qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - &
    1813                &          (sin36*b(ia+i)+sin72*b(ib+i))
    1814           c(je+j) = -(qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - &
    1815                &          (sin36*b(ia+i)+sin72*b(ib+i))
    1816           c(jc+j) = (qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - &
    1817                &          (sin72*b(ia+i)-sin36*b(ib+i))
    1818           c(jd+j) = -(qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) - &
    1819                &          (sin72*b(ia+i)-sin36*b(ib+i))
    1820           i = i + inc3
    1821           j = j + inc4
    1822        END DO
    1823        ibase = ibase + inc1
    1824        jbase = jbase + inc2
    1825     END DO
    1826 
    1827     GO TO 170
    1828 120 CONTINUE
    1829     qqrt5 = 2.0_wp*qrt5
    1830     ssin36 = 2.0_wp*sin36
    1831     ssin72 = 2.0_wp*sin72
    1832     DO l = 1, la
    1833        i = ibase
    1834        j = jbase
    1835 !DIR$ IVDEP
    1836 !CDIR NODEP
    1837 !OCL NOVREC
    1838        DO ijk = 1, lot
    1839           c(ja+j) = 2.0_wp*(a(ia+i)+(a(ib+i)+a(ic+i)))
    1840           c(jb+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ &
    1841                &          i))) - (ssin72*b(ib+i)+ssin36*b(ic+i))
    1842           c(jc+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ &
    1843                &          i))) - (ssin36*b(ib+i)-ssin72*b(ic+i))
    1844           c(jd+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qqrt5*(a(ib+i)-a(ic+ &
    1845                &          i))) + (ssin36*b(ib+i)-ssin72*b(ic+i))
    1846           c(je+j) = (2.0_wp*(a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qqrt5*(a(ib+i)-a(ic+ &
    1847                &          i))) + (ssin72*b(ib+i)+ssin36*b(ic+i))
    1848           i = i + inc3
    1849           j = j + inc4
    1850        END DO
    1851        ibase = ibase + inc1
    1852        jbase = jbase + inc2
    1853     END DO
    1854     GO TO 170
    1855 
    1856     ! Coding for factor 6
    1857 
    1858 130 CONTINUE
    1859     ia = 1
    1860     ib = ia + (2*m-la)*inc1
    1861     ic = ib + 2*m*inc1
    1862     id = ic + 2*m*inc1
    1863     ie = ic
    1864     if = ib
    1865     ja = 1
    1866     jb = ja + jink
    1867     jc = jb + jink
    1868     jd = jc + jink
    1869     je = jd + jink
    1870     jf = je + jink
    1871 
    1872     IF (la==m) GO TO 150
    1873 
    1874     DO l = 1, la
    1875        i = ibase
    1876        j = jbase
    1877 !DIR$ IVDEP
    1878 !CDIR NODEP
    1879 !OCL NOVREC
    1880        DO ijk = 1, lot
    1881           c(ja+j) = (a(ia+i)+a(id+i)) + (a(ib+i)+a(ic+i))
    1882           c(jd+j) = (a(ia+i)-a(id+i)) - (a(ib+i)-a(ic+i))
    1883           c(jb+j) = ((a(ia+i)-a(id+i))+0.5_wp*(a(ib+i)-a(ic+i))) - (sin60*(b(ib+ &
    1884                &          i)+b(ic+i)))
    1885           c(jf+j) = ((a(ia+i)-a(id+i))+0.5_wp*(a(ib+i)-a(ic+i))) + (sin60*(b(ib+ &
    1886                &          i)+b(ic+i)))
    1887           c(jc+j) = ((a(ia+i)+a(id+i))-0.5_wp*(a(ib+i)+a(ic+i))) - (sin60*(b(ib+ &
    1888                &          i)-b(ic+i)))
    1889           c(je+j) = ((a(ia+i)+a(id+i))-0.5_wp*(a(ib+i)+a(ic+i))) + (sin60*(b(ib+ &
    1890                &          i)-b(ic+i)))
    1891           i = i + inc3
    1892           j = j + inc4
    1893        END DO
    1894        ibase = ibase + inc1
    1895        jbase = jbase + inc2
    1896     END DO
    1897     ia = ia + iink
    1898     iink = 2*iink
    1899     ib = ib + iink
    1900     ic = ic + iink
    1901     id = id - iink
    1902     ie = ie - iink
    1903     if = if - iink
    1904     jbase = jbase + jump
    1905     jump = 2*jump + jink
    1906     IF (ic==id) GO TO 140
    1907     DO k = la, kstop, la
    1908        kb = k + k
    1909        kc = kb + kb
    1910        kd = kc + kb
    1911        ke = kd + kb
    1912        kf = ke + kb
    1913        c1 = trigs(kb+1)
    1914        s1 = trigs(kb+2)
    1915        c2 = trigs(kc+1)
    1916        s2 = trigs(kc+2)
    1917        c3 = trigs(kd+1)
    1918        s3 = trigs(kd+2)
    1919        c4 = trigs(ke+1)
    1920        s4 = trigs(ke+2)
    1921        c5 = trigs(kf+1)
    1922        s5 = trigs(kf+2)
    1923        ibase = 0
    1924        DO l = 1, la
    1925           i = ibase
    1926           j = jbase
    1927 !DIR$ IVDEP
    1928 !CDIR NODEP
    1929 !OCL NOVREC
    1930           DO ijk = 1, lot
    1931 
    1932              a11(ijk) = (a(ie+i)+a(ib+i)) + (a(ic+i)+a(if+i))
    1933              a20(ijk) = (a(ia+i)+a(id+i)) - 0.5_wp*a11(ijk)
    1934              a21(ijk) = sin60*((a(ie+i)+a(ib+i))-(a(ic+i)+a(if+i)))
    1935              b11(ijk) = (b(ib+i)-b(ie+i)) + (b(ic+i)-b(if+i))
    1936              b20(ijk) = (b(ia+i)-b(id+i)) - 0.5_wp*b11(ijk)
    1937              b21(ijk) = sin60*((b(ib+i)-b(ie+i))-(b(ic+i)-b(if+i)))
    1938 
    1939              c(ja+j) = (a(ia+i)+a(id+i)) + a11(ijk)
    1940              d(ja+j) = (b(ia+i)-b(id+i)) + b11(ijk)
    1941              c(jc+j) = c2*(a20(ijk)-b21(ijk)) - s2*(b20(ijk)+a21(ijk))
    1942              d(jc+j) = s2*(a20(ijk)-b21(ijk)) + c2*(b20(ijk)+a21(ijk))
    1943              c(je+j) = c4*(a20(ijk)+b21(ijk)) - s4*(b20(ijk)-a21(ijk))
    1944              d(je+j) = s4*(a20(ijk)+b21(ijk)) + c4*(b20(ijk)-a21(ijk))
    1945 
    1946              a11(ijk) = (a(ie+i)-a(ib+i)) + (a(ic+i)-a(if+i))
    1947              b11(ijk) = (b(ie+i)+b(ib+i)) - (b(ic+i)+b(if+i))
    1948              a20(ijk) = (a(ia+i)-a(id+i)) - 0.5_wp*a11(ijk)
    1949              a21(ijk) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
    1950              b20(ijk) = (b(ia+i)+b(id+i)) + 0.5_wp*b11(ijk)
    1951              b21(ijk) = sin60*((b(ie+i)+b(ib+i))+(b(ic+i)+b(if+i)))
    1952 
    1953              c(jd+j) = c3*((a(ia+i)-a(id+i))+a11(ijk)) - s3*((b(ia+i)+b(id+ &
    1954                   &            i))-b11(ijk))
    1955              d(jd+j) = s3*((a(ia+i)-a(id+i))+a11(ijk)) + c3*((b(ia+i)+b(id+ &
    1956                   &            i))-b11(ijk))
    1957              c(jb+j) = c1*(a20(ijk)-b21(ijk)) - s1*(b20(ijk)-a21(ijk))
    1958              d(jb+j) = s1*(a20(ijk)-b21(ijk)) + c1*(b20(ijk)-a21(ijk))
    1959              c(jf+j) = c5*(a20(ijk)+b21(ijk)) - s5*(b20(ijk)+a21(ijk))
    1960              d(jf+j) = s5*(a20(ijk)+b21(ijk)) + c5*(b20(ijk)+a21(ijk))
    1961 
    1962              i = i + inc3
    1963              j = j + inc4
    1964           END DO
    1965           ibase = ibase + inc1
    1966           jbase = jbase + inc2
    1967        END DO
    1968        ia = ia + iink
    1969        ib = ib + iink
    1970        ic = ic + iink
    1971        id = id - iink
    1972        ie = ie - iink
    1973        if = if - iink
    1974        jbase = jbase + jump
    1975     END DO
    1976     IF (ic>id) GO TO 170
    1977 140 CONTINUE
    1978     ibase = 0
    1979     DO l = 1, la
    1980        i = ibase
    1981        j = jbase
    1982 !DIR$ IVDEP
    1983 !CDIR NODEP
    1984 !OCL NOVREC
    1985        DO ijk = 1, lot
    1986           c(ja+j) = a(ib+i) + (a(ia+i)+a(ic+i))
    1987           c(jd+j) = b(ib+i) - (b(ia+i)+b(ic+i))
    1988           c(jb+j) = (sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i))
    1989           c(jf+j) = -(sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i))
    1990           c(jc+j) = sin60*(b(ic+i)-b(ia+i)) + (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i))
    1991           c(je+j) = sin60*(b(ic+i)-b(ia+i)) - (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i))
    1992           i = i + inc3
    1993           j = j + inc4
    1994        END DO
    1995        ibase = ibase + inc1
    1996        jbase = jbase + inc2
    1997     END DO
    1998 
    1999     GO TO 170
    2000 150 CONTINUE
    2001     ssin60 = 2.0_wp*sin60
    2002     DO l = 1, la
    2003        i = ibase
    2004        j = jbase
    2005 !DIR$ IVDEP
    2006 !CDIR NODEP
    2007 !OCL NOVREC
    2008        DO ijk = 1, lot
    2009           c(ja+j) = (2.0_wp*(a(ia+i)+a(id+i))) + (2.0_wp*(a(ib+i)+a(ic+i)))
    2010           c(jd+j) = (2.0_wp*(a(ia+i)-a(id+i))) - (2.0_wp*(a(ib+i)-a(ic+i)))
    2011           c(jb+j) = (2.0_wp*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) - (ssin60*(b(ib+ &
    2012                &          i)+b(ic+i)))
    2013           c(jf+j) = (2.0_wp*(a(ia+i)-a(id+i))+(a(ib+i)-a(ic+i))) + (ssin60*(b(ib+ &
    2014                &          i)+b(ic+i)))
    2015           c(jc+j) = (2.0_wp*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) - (ssin60*(b(ib+ &
    2016                &          i)-b(ic+i)))
    2017           c(je+j) = (2.0_wp*(a(ia+i)+a(id+i))-(a(ib+i)+a(ic+i))) + (ssin60*(b(ib+ &
    2018                &          i)-b(ic+i)))
    2019           i = i + inc3
    2020           j = j + inc4
    2021        END DO
    2022        ibase = ibase + inc1
    2023        jbase = jbase + inc2
    2024     END DO
    2025     GO TO 170
    2026 
    2027     ! Coding for factor 8
    2028 
    2029 160 CONTINUE
    2030     ibad = 3
    2031     IF (la/=m) GO TO 180
    2032     ia = 1
    2033     ib = ia + la*inc1
    2034     ic = ib + 2*la*inc1
    2035     id = ic + 2*la*inc1
    2036     ie = id + 2*la*inc1
    2037     ja = 1
    2038     jb = ja + jink
    2039     jc = jb + jink
    2040     jd = jc + jink
    2041     je = jd + jink
    2042     jf = je + jink
    2043     jg = jf + jink
    2044     jh = jg + jink
    2045     ssin45 = SQRT(2.0_wp)
    2046 
    2047     DO l = 1, la
    2048        i = ibase
    2049        j = jbase
    2050 !DIR$ IVDEP
    2051 !CDIR NODEP
    2052 !OCL NOVREC
    2053        DO ijk = 1, lot
    2054           c(ja+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))+(a(ib+i)+a(id+i)))
    2055           c(je+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))-(a(ib+i)+a(id+i)))
    2056           c(jc+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))-(b(ib+i)-b(id+i)))
    2057           c(jg+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))+(b(ib+i)-b(id+i)))
    2058           c(jb+j) = 2.0_wp*((a(ia+i)-a(ie+i))-b(ic+i)) + ssin45*((a(ib+i)-a(id+ &
    2059                &              i))-(b(ib+i)+b(id+i)))
    2060           c(jf+j) = 2.0_wp*((a(ia+i)-a(ie+i))-b(ic+i)) - ssin45*((a(ib+i)-a(id+ &
    2061                &              i))-(b(ib+i)+b(id+i)))
    2062           c(jd+j) = 2.0_wp*((a(ia+i)-a(ie+i))+b(ic+i)) - ssin45*((a(ib+i)-a(id+ &
    2063                &              i))+(b(ib+i)+b(id+i)))
    2064           c(jh+j) = 2.0_wp*((a(ia+i)-a(ie+i))+b(ic+i)) + ssin45*((a(ib+i)-a(id+ &
    2065                &              i))+(b(ib+i)+b(id+i)))
    2066           i = i + inc3
    2067           j = j + inc4
    2068        END DO
    2069        ibase = ibase + inc1
    2070        jbase = jbase + inc2
    2071     END DO
    2072 
    2073     ! Return
    2074 
    2075 170 CONTINUE
    2076     ibad = 0
    2077 180 CONTINUE
    2078     ierr = ibad
    2079     RETURN
    2080   END SUBROUTINE rpassm
     1387    IF ( igo == 7 )  igo = 6
     1388    IF ( igo < 1  .OR.  igo > 6 )  THEN
     1389       ierr = 2
     1390       RETURN
     1391    ENDIF
     1392
     1393
     1394    SELECT CASE ( igo )
     1395!
     1396!--    Coding for factor 2
     1397       CASE ( 1 )
     1398
     1399          ia = 1
     1400          ib = ia + (2*m-la) * inc1
     1401          ja = 1
     1402          jb = ja + jink
     1403
     1404          IF ( la /= m )  THEN
     1405
     1406             DO  l = 1, la
     1407                i = ibase
     1408                j = jbase
     1409                !DIR$ IVDEP
     1410                !CDIR NODEP
     1411                !OCL NOVREC
     1412                DO  ijk = 1, lot
     1413                   c(ja+j) = a(ia+i) + a(ib+i)
     1414                   c(jb+j) = a(ia+i) - a(ib+i)
     1415                   i = i + inc3
     1416                   j = j + inc4
     1417                ENDDO
     1418                ibase = ibase + inc1
     1419                jbase = jbase + inc2
     1420             ENDDO
     1421
     1422             ia = ia + iink
     1423             iink = 2 * iink
     1424             ib = ib - iink
     1425             ibase = 0
     1426             jbase = jbase + jump
     1427             jump = 2 * jump + jink
     1428
     1429             IF ( ia /= ib )  THEN
     1430
     1431                DO  k = la, kstop, la
     1432                   kb = k + k
     1433                   c1 = trigs(kb+1)
     1434                   s1 = trigs(kb+2)
     1435                   ibase = 0
     1436                   DO  l = 1, la
     1437                      i = ibase
     1438                      j = jbase
     1439                      !DIR$ IVDEP
     1440                      !CDIR NODEP
     1441                      !OCL NOVREC
     1442                      DO  ijk = 1, lot
     1443                         c(ja+j) = a(ia+i) + a(ib+i)
     1444                         d(ja+j) = b(ia+i) - b(ib+i)
     1445                         c(jb+j) = c1*(a(ia+i)-a(ib+i)) - s1*(b(ia+i)+b(ib+i))
     1446                         d(jb+j) = s1*(a(ia+i)-a(ib+i)) + c1*(b(ia+i)+b(ib+i))
     1447                         i = i + inc3
     1448                         j = j + inc4
     1449                      ENDDO
     1450                      ibase = ibase + inc1
     1451                      jbase = jbase + inc2
     1452                   ENDDO
     1453                   ia = ia + iink
     1454                   ib = ib - iink
     1455                   jbase = jbase + jump
     1456                ENDDO
     1457
     1458                IF ( ia > ib )  RETURN
     1459
     1460             ENDIF
     1461
     1462             ibase = 0
     1463             DO  l = 1, la
     1464                i = ibase
     1465                j = jbase
     1466                !DIR$ IVDEP
     1467                !CDIR NODEP
     1468                !OCL NOVREC
     1469                DO  ijk = 1, lot
     1470                   c(ja+j) = a(ia+i)
     1471                   c(jb+j) = -b(ia+i)
     1472                   i = i + inc3
     1473                   j = j + inc4
     1474                ENDDO
     1475                ibase = ibase + inc1
     1476                jbase = jbase + inc2
     1477             ENDDO
     1478
     1479          ELSE
     1480
     1481             DO  l = 1, la
     1482                i = ibase
     1483                j = jbase
     1484                !DIR$ IVDEP
     1485                !CDIR NODEP
     1486                !OCL NOVREC
     1487                DO  ijk = 1, lot
     1488                   c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i))
     1489                   c(jb+j) = 2.0_wp*(a(ia+i)-a(ib+i))
     1490                   i = i + inc3
     1491                   j = j + inc4
     1492                ENDDO
     1493                ibase = ibase + inc1
     1494                jbase = jbase + inc2
     1495             ENDDO
     1496
     1497          ENDIF
     1498
     1499!
     1500!--    Coding for factor 3
     1501       CASE ( 2 )
     1502
     1503          ia = 1
     1504          ib = ia + (2*m-la) * inc1
     1505          ic = ib
     1506          ja = 1
     1507          jb = ja + jink
     1508          jc = jb + jink
     1509
     1510          IF ( la /= m )  THEN
     1511
     1512             DO  l = 1, la
     1513                i = ibase
     1514                j = jbase
     1515                !DIR$ IVDEP
     1516                !CDIR NODEP
     1517                !OCL NOVREC
     1518                DO  ijk = 1, lot
     1519                   c(ja+j) = a(ia+i) + a(ib+i)
     1520                   c(jb+j) = (a(ia+i)-0.5_wp*a(ib+i)) - (sin60*(b(ib+i)))
     1521                   c(jc+j) = (a(ia+i)-0.5_wp*a(ib+i)) + (sin60*(b(ib+i)))
     1522                   i = i + inc3
     1523                   j = j + inc4
     1524                ENDDO
     1525                ibase = ibase + inc1
     1526                jbase = jbase + inc2
     1527             ENDDO
     1528             ia = ia + iink
     1529             iink = 2 * iink
     1530             ib = ib + iink
     1531             ic = ic - iink
     1532             jbase = jbase + jump
     1533             jump = 2 * jump + jink
     1534
     1535             IF ( ia /= ic )  THEN
     1536
     1537                DO  k = la, kstop, la
     1538                   kb = k + k
     1539                   kc = kb + kb
     1540                   c1 = trigs(kb+1)
     1541                   s1 = trigs(kb+2)
     1542                   c2 = trigs(kc+1)
     1543                   s2 = trigs(kc+2)
     1544                   ibase = 0
     1545                   DO  l = 1, la
     1546                      i = ibase
     1547                      j = jbase
     1548                      !DIR$ IVDEP
     1549                      !CDIR NODEP
     1550                      !OCL NOVREC
     1551                      DO  ijk = 1, lot
     1552                         c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
     1553                         d(ja+j) = b(ia+i) + (b(ib+i)-b(ic+i))
     1554                         c(jb+j) = c1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ b(ic+i)))) &
     1555                                   - s1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- a(ic+i))))
     1556                         d(jb+j) = s1*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))-(sin60*(b(ib+i)+ b(ic+i)))) &
     1557                                   + c1*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))+(sin60*(a(ib+i)- a(ic+i))))
     1558                         c(jc+j) = c2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ b(ic+i)))) &
     1559                                   - s2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- a(ic+i))))
     1560                         d(jc+j) = s2*((a(ia+i)-0.5_wp*(a(ib+i)+a(ic+i)))+(sin60*(b(ib+i)+ b(ic+i)))) &
     1561                                   + c2*((b(ia+i)-0.5_wp*(b(ib+i)-b(ic+i)))-(sin60*(a(ib+i)- a(ic+i))))
     1562                         i = i + inc3
     1563                         j = j + inc4
     1564                      ENDDO
     1565                      ibase = ibase + inc1
     1566                      jbase = jbase + inc2
     1567                   ENDDO
     1568                   ia = ia + iink
     1569                   ib = ib + iink
     1570                   ic = ic - iink
     1571                   jbase = jbase + jump
     1572                ENDDO
     1573
     1574                IF ( ia > ic )  RETURN
     1575
     1576             ENDIF
     1577
     1578             ibase = 0
     1579             DO  l = 1, la
     1580                i = ibase
     1581                j = jbase
     1582                !DIR$ IVDEP
     1583                !CDIR NODEP
     1584                !OCL NOVREC
     1585                DO  ijk = 1, lot
     1586                   c(ja+j) = a(ia+i) + a(ib+i)
     1587                   c(jb+j) = (0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
     1588                   c(jc+j) = -(0.5_wp*a(ia+i)-a(ib+i)) - (sin60*b(ia+i))
     1589                   i = i + inc3
     1590                   j = j + inc4
     1591                ENDDO
     1592                ibase = ibase + inc1
     1593                jbase = jbase + inc2
     1594             ENDDO
     1595
     1596          ELSE
     1597
     1598             ssin60 = 2.0_wp * sin60
     1599             DO  l = 1, la
     1600                i = ibase
     1601                j = jbase
     1602                !DIR$ IVDEP
     1603                !CDIR NODEP
     1604                !OCL NOVREC
     1605                DO  ijk = 1, lot
     1606                   c(ja+j) = 2.0_wp*(a(ia+i)+a(ib+i))
     1607                   c(jb+j) = (2.0_wp*a(ia+i)-a(ib+i)) - (ssin60*b(ib+i))
     1608                   c(jc+j) = (2.0_wp*a(ia+i)-a(ib+i)) + (ssin60*b(ib+i))
     1609                   i = i + inc3
     1610                   j = j + inc4
     1611                ENDDO
     1612                ibase = ibase + inc1
     1613                jbase = jbase + inc2
     1614             ENDDO
     1615
     1616          ENDIF
     1617
     1618!
     1619!--    Coding for factor 4
     1620       CASE ( 3 )
     1621
     1622          ia = 1
     1623          ib = ia + (2*m-la) * inc1
     1624          ic = ib + 2*m*inc1
     1625          id = ib
     1626          ja = 1
     1627          jb = ja + jink
     1628          jc = jb + jink
     1629          jd = jc + jink
     1630
     1631          IF ( la /= m )  THEN
     1632
     1633             DO  l = 1, la
     1634                i = ibase
     1635                j = jbase
     1636                !DIR$ IVDEP
     1637                !CDIR NODEP
     1638                !OCL NOVREC
     1639                DO  ijk = 1, lot
     1640                   c(ja+j) = (a(ia+i)+a(ic+i)) + a(ib+i)
     1641                   c(jb+j) = (a(ia+i)-a(ic+i)) - b(ib+i)
     1642                   c(jc+j) = (a(ia+i)+a(ic+i)) - a(ib+i)
     1643                   c(jd+j) = (a(ia+i)-a(ic+i)) + b(ib+i)
     1644                   i = i + inc3
     1645                   j = j + inc4
     1646               ENDDO
     1647                ibase = ibase + inc1
     1648                jbase = jbase + inc2
     1649             ENDDO
     1650             ia = ia + iink
     1651             iink = 2 * iink
     1652             ib = ib + iink
     1653             ic = ic - iink
     1654             id = id - iink
     1655             jbase = jbase + jump
     1656             jump = 2 * jump + jink
     1657
     1658             IF ( ib /= ic )  THEN
     1659
     1660                DO  k = la, kstop, la
     1661                   kb = k + k
     1662                   kc = kb + kb
     1663                   kd = kc + kb
     1664                   c1 = trigs(kb+1)
     1665                   s1 = trigs(kb+2)
     1666                   c2 = trigs(kc+1)
     1667                   s2 = trigs(kc+2)
     1668                   c3 = trigs(kd+1)
     1669                   s3 = trigs(kd+2)
     1670                   ibase = 0
     1671                   DO  l = 1, la
     1672                      i = ibase
     1673                      j = jbase
     1674                      !DIR$ IVDEP
     1675                      !CDIR NODEP
     1676                      !OCL NOVREC
     1677                      DO  ijk = 1, lot
     1678                         c(ja+j) = (a(ia+i)+a(ic+i)) + (a(ib+i)+a(id+i))
     1679                         d(ja+j) = (b(ia+i)-b(ic+i)) + (b(ib+i)-b(id+i))
     1680                         c(jc+j) = c2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) - s2*((b(ia+i)-b(ic+i))&
     1681                                   -(b(ib+i)-b(id+i)))
     1682                         d(jc+j) = s2*((a(ia+i)+a(ic+i))-(a(ib+i)+a(id+i))) + c2*((b(ia+i)-b(ic+i))&
     1683                                   -(b(ib+i)-b(id+i)))
     1684                         c(jb+j) = c1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) - s1*((b(ia+i)+b(ic+i))&
     1685                                   +(a(ib+i)-a(id+i)))
     1686                         d(jb+j) = s1*((a(ia+i)-a(ic+i))-(b(ib+i)+b(id+i))) + c1*((b(ia+i)+b(ic+i))&
     1687                                   +(a(ib+i)-a(id+i)))
     1688                         c(jd+j) = c3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) - s3*((b(ia+i)+b(ic+i))&
     1689                                   -(a(ib+i)-a(id+i)))
     1690                         d(jd+j) = s3*((a(ia+i)-a(ic+i))+(b(ib+i)+b(id+i))) + c3*((b(ia+i)+b(ic+i))&
     1691                                   -(a(ib+i)-a(id+i)))
     1692                         i = i + inc3
     1693                         j = j + inc4
     1694                      ENDDO
     1695                      ibase = ibase + inc1
     1696                      jbase = jbase + inc2
     1697                   ENDDO
     1698                   ia = ia + iink
     1699                   ib = ib + iink
     1700                   ic = ic - iink
     1701                   id = id - iink
     1702                   jbase = jbase + jump
     1703                ENDDO
     1704
     1705                IF ( ib > ic )  RETURN
     1706
     1707             ENDIF
     1708
     1709             ibase = 0
     1710             sin45 = SQRT( 0.5_wp )
     1711             DO  l = 1, la
     1712                i = ibase
     1713                j = jbase
     1714                !DIR$ IVDEP
     1715                !CDIR NODEP
     1716                !OCL NOVREC
     1717                DO  ijk = 1, lot
     1718                   c(ja+j) = a(ia+i) + a(ib+i)
     1719                   c(jb+j) = sin45*((a(ia+i)-a(ib+i))-(b(ia+i)+b(ib+i)))
     1720                   c(jc+j) = b(ib+i) - b(ia+i)
     1721                   c(jd+j) = -sin45*((a(ia+i)-a(ib+i))+(b(ia+i)+b(ib+i)))
     1722                   i = i + inc3
     1723                   j = j + inc4
     1724                ENDDO
     1725                ibase = ibase + inc1
     1726                jbase = jbase + inc2
     1727             ENDDO
     1728
     1729          ELSE
     1730
     1731             DO  l = 1, la
     1732                i = ibase
     1733                j = jbase
     1734                !DIR$ IVDEP
     1735                !CDIR NODEP
     1736                !OCL NOVREC
     1737                DO  ijk = 1, lot
     1738                   c(ja+j) = 2.0_wp*((a(ia+i)+a(ic+i))+a(ib+i))
     1739                   c(jb+j) = 2.0_wp*((a(ia+i)-a(ic+i))-b(ib+i))
     1740                   c(jc+j) = 2.0_wp*((a(ia+i)+a(ic+i))-a(ib+i))
     1741                   c(jd+j) = 2.0_wp*((a(ia+i)-a(ic+i))+b(ib+i))
     1742                   i = i + inc3
     1743                   j = j + inc4
     1744                ENDDO
     1745                ibase = ibase + inc1
     1746                jbase = jbase + inc2
     1747             ENDDO
     1748
     1749          ENDIF
     1750
     1751!
     1752!--    Coding for factor 5
     1753       CASE ( 4 )
     1754
     1755          ia = 1
     1756          ib = ia + (2*m-la) * inc1
     1757          ic = ib + 2*m*inc1
     1758          id = ic
     1759          ie = ib
     1760          ja = 1
     1761          jb = ja + jink
     1762          jc = jb + jink
     1763          jd = jc + jink
     1764          je = jd + jink
     1765
     1766          IF ( la /= m )  THEN
     1767
     1768             DO  l = 1, la
     1769                i = ibase
     1770                j = jbase
     1771                !DIR$ IVDEP
     1772                !CDIR NODEP
     1773                !OCL NOVREC
     1774                DO  ijk = 1, lot
     1775                   c(ja+j) = a(ia+i) + (a(ib+i)+a(ic+i))
     1776                   c(jb+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) -        &
     1777                             (sin72*b(ib+i)+sin36*b(ic+i))
     1778                   c(jc+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) -        &
     1779                             (sin36*b(ib+i)-sin72*b(ic+i))
     1780                   c(jd+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))-qrt5*(a(ib+i)-a(ic+i))) +        &
     1781                             (sin36*b(ib+i)-sin72*b(ic+i))
     1782                   c(je+j) = ((a(ia+i)-0.25_wp*(a(ib+i)+a(ic+i)))+qrt5*(a(ib+i)-a(ic+i))) +        &
     1783                             (sin72*b(ib+i)+sin36*b(ic+i))
     1784                   i = i + inc3
     1785                   j = j + inc4
     1786                ENDDO
     1787                ibase = ibase + inc1
     1788                jbase = jbase + inc2
     1789             ENDDO
     1790             ia = ia + iink
     1791             iink = 2 * iink
     1792             ib = ib + iink
     1793             ic = ic + iink
     1794             id = id - iink
     1795             ie = ie - iink
     1796             jbase = jbase + jump
     1797             jump = 2 * jump + jink
     1798
     1799             IF ( ib /= id )  THEN
     1800
     1801                DO  k = la, kstop, la
     1802                   kb = k + k
     1803                   kc = kb + kb
     1804                   kd = kc + kb
     1805                   ke = kd + kb
     1806                   c1 = trigs(kb+1)
     1807                   s1 = trigs(kb+2)
     1808                   c2 = trigs(kc+1)
     1809                   s2 = trigs(kc+2)
     1810                   c3 = trigs(kd+1)
     1811                   s3 = trigs(kd+2)
     1812                   c4 = trigs(ke+1)
     1813                   s4 = trigs(ke+2)
     1814                   ibase = 0
     1815                   DO  l = 1, la
     1816                      i = ibase
     1817                      j = jbase
     1818                      !DIR$ IVDEP
     1819                      !CDIR NODEP
     1820                      !OCL NOVREC
     1821                      DO  ijk = 1, lot
     1822                         a10(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) +      &
     1823                                    qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
     1824                         a20(ijk) = (a(ia+i)-0.25_wp*((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))) -      &
     1825                                    qrt5*((a(ib+i)+a(ie+i))-(a(ic+i)+a(id+i)))
     1826                         b10(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) +      &
     1827                                    qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
     1828                         b20(ijk) = (b(ia+i)-0.25_wp*((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))) -      &
     1829                                    qrt5*((b(ib+i)-b(ie+i))-(b(ic+i)-b(id+i)))
     1830                         a11(ijk) = sin72*(b(ib+i)+b(ie+i)) + sin36*(b(ic+i)+b(id+i))
     1831                         a21(ijk) = sin36*(b(ib+i)+b(ie+i)) - sin72*(b(ic+i)+b(id+i))
     1832                         b11(ijk) = sin72*(a(ib+i)-a(ie+i)) + sin36*(a(ic+i)-a(id+i))
     1833                         b21(ijk) = sin36*(a(ib+i)-a(ie+i)) - sin72*(a(ic+i)-a(id+i))
     1834                         c(ja+j) = a(ia+i) + ((a(ib+i)+a(ie+i))+(a(ic+i)+a(id+i)))
     1835                         d(ja+j) = b(ia+i) + ((b(ib+i)-b(ie+i))+(b(ic+i)-b(id+i)))
     1836                         c(jb+j) = c1*(a10(ijk)-a11(ijk)) - s1*(b10(ijk)+b11(ijk))
     1837                         d(jb+j) = s1*(a10(ijk)-a11(ijk)) + c1*(b10(ijk)+b11(ijk))
     1838                         c(je+j) = c4*(a10(ijk)+a11(ijk)) - s4*(b10(ijk)-b11(ijk))
     1839                         d(je+j) = s4*(a10(ijk)+a11(ijk)) + c4*(b10(ijk)-b11(ijk))
     1840                         c(jc+j) = c2*(a20(ijk)-a21(ijk)) - s2*(b20(ijk)+b21(ijk))
     1841                         d(jc+j) = s2*(a20(ijk)-a21(ijk)) + c2*(b20(ijk)+b21(ijk))
     1842                         c(jd+j) = c3*(a20(ijk)+a21(ijk)) - s3*(b20(ijk)-b21(ijk))
     1843                         d(jd+j) = s3*(a20(ijk)+a21(ijk)) + c3*(b20(ijk)-b21(ijk))
     1844
     1845                         i = i + inc3
     1846                         j = j + inc4
     1847                      ENDDO
     1848                      ibase = ibase + inc1
     1849                      jbase = jbase + inc2
     1850                   ENDDO
     1851                   ia = ia + iink
     1852                   ib = ib + iink
     1853                   ic = ic + iink
     1854                   id = id - iink
     1855                   ie = ie - iink
     1856                   jbase = jbase + jump
     1857                ENDDO
     1858
     1859                IF ( ib > id )  RETURN
     1860
     1861             ENDIF
     1862
     1863             ibase = 0
     1864             DO  l = 1, la
     1865                i = ibase
     1866                j = jbase
     1867                !DIR$ IVDEP
     1868                !CDIR NODEP
     1869                !OCL NOVREC
     1870                DO  ijk = 1, lot
     1871                   c(ja+j) = (a(ia+i)+a(ib+i)) + a(ic+i)
     1872                   c(jb+j) = (qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) -        &
     1873                             (sin36*b(ia+i)+sin72*b(ib+i))
     1874                   c(je+j) = -(qrt5*(a(ia+i)-a(ib+i))+(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) -       &
     1875                             (sin36*b(ia+i)+sin72*b(ib+i))
     1876                   c(jc+j) = (qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) -        &
     1877                             (sin72*b(ia+i)-sin36*b(ib+i))
     1878                   c(jd+j) = -(qrt5*(a(ia+i)-a(ib+i))-(0.25_wp*(a(ia+i)+a(ib+i))-a(ic+i))) -       &
     1879                             (sin72*b(ia+i)-sin36*b(ib+i))
     1880                   i = i + inc3
     1881                   j = j + inc4
     1882                ENDDO
     1883                ibase = ibase + inc1
     1884                jbase = jbase + inc2
     1885             ENDDO
     1886
     1887          ELSE
     1888
     1889             qqrt5  = 2.0_wp * qrt5
     1890             ssin36 = 2.0_wp * sin36
     1891             ssin72 = 2.0_wp * sin72
     1892             DO  l = 1, la
     1893                i = ibase
     1894                j = jbase
     1895                !DIR$ IVDEP
     1896                !CDIR NODEP
     1897                !OCL NOVREC
     1898                DO  ijk = 1, lot
     1899                   c(ja+j) = 2.0_wp*(a(ia+i)+(a(ib+i)+a(ic+i)))
     1900                   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)))  &
     1901                             - (ssin72*b(ib+i)+ssin36*b(ic+i))
     1902                   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)))  &
     1903                             - (ssin36*b(ib+i)-ssin72*b(ic+i))
     1904                   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)))  &
     1905                             + (ssin36*b(ib+i)-ssin72*b(ic+i))
     1906                   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)))  &
     1907                             + (ssin72*b(ib+i)+ssin36*b(ic+i))
     1908                   i = i + inc3
     1909                   j = j + inc4
     1910                ENDDO
     1911                ibase = ibase + inc1
     1912                jbase = jbase + inc2
     1913             ENDDO
     1914
     1915          ENDIF
     1916
     1917!
     1918!--    Coding for factor 6
     1919       CASE ( 5 )
     1920
     1921          ia = 1
     1922          ib = ia + (2*m-la) * inc1
     1923          ic = ib + 2*m*inc1
     1924          id = ic + 2*m*inc1
     1925          ie = ic
     1926          if = ib
     1927          ja = 1
     1928          jb = ja + jink
     1929          jc = jb + jink
     1930          jd = jc + jink
     1931          je = jd + jink
     1932          jf = je + jink
     1933
     1934          IF ( la /= m )  THEN
     1935
     1936             DO  l = 1, la
     1937                i = ibase
     1938                j = jbase
     1939                !DIR$ IVDEP
     1940                !CDIR NODEP
     1941                !OCL NOVREC
     1942                DO  ijk = 1, lot
     1943                   c(ja+j) = (a(ia+i)+a(id+i)) + (a(ib+i)+a(ic+i))
     1944                   c(jd+j) = (a(ia+i)-a(id+i)) - (a(ib+i)-a(ic+i))
     1945                   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)))
     1946                   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)))
     1947                   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)))
     1948                   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)))
     1949                   i = i + inc3
     1950                   j = j + inc4
     1951                ENDDO
     1952                ibase = ibase + inc1
     1953                jbase = jbase + inc2
     1954             ENDDO
     1955             ia = ia + iink
     1956             iink = 2 * iink
     1957             ib = ib + iink
     1958             ic = ic + iink
     1959             id = id - iink
     1960             ie = ie - iink
     1961             if = if - iink
     1962             jbase = jbase + jump
     1963             jump = 2 * jump + jink
     1964
     1965             IF ( ic /= id )  THEN
     1966
     1967                DO  k = la, kstop, la
     1968                   kb = k + k
     1969                   kc = kb + kb
     1970                   kd = kc + kb
     1971                   ke = kd + kb
     1972                   kf = ke + kb
     1973                   c1 = trigs(kb+1)
     1974                   s1 = trigs(kb+2)
     1975                   c2 = trigs(kc+1)
     1976                   s2 = trigs(kc+2)
     1977                   c3 = trigs(kd+1)
     1978                   s3 = trigs(kd+2)
     1979                   c4 = trigs(ke+1)
     1980                   s4 = trigs(ke+2)
     1981                   c5 = trigs(kf+1)
     1982                   s5 = trigs(kf+2)
     1983                   ibase = 0
     1984                   DO  l = 1, la
     1985                      i = ibase
     1986                      j = jbase
     1987                      !DIR$ IVDEP
     1988                      !CDIR NODEP
     1989                      !OCL NOVREC
     1990                      DO  ijk = 1, lot
     1991                         a11(ijk) = (a(ie+i)+a(ib+i)) + (a(ic+i)+a(if+i))
     1992                         a20(ijk) = (a(ia+i)+a(id+i)) - 0.5_wp*a11(ijk)
     1993                         a21(ijk) = sin60*((a(ie+i)+a(ib+i))-(a(ic+i)+a(if+i)))
     1994                         b11(ijk) = (b(ib+i)-b(ie+i)) + (b(ic+i)-b(if+i))
     1995                         b20(ijk) = (b(ia+i)-b(id+i)) - 0.5_wp*b11(ijk)
     1996                         b21(ijk) = sin60*((b(ib+i)-b(ie+i))-(b(ic+i)-b(if+i)))
     1997
     1998                         c(ja+j) = (a(ia+i)+a(id+i)) + a11(ijk)
     1999                         d(ja+j) = (b(ia+i)-b(id+i)) + b11(ijk)
     2000                         c(jc+j) = c2*(a20(ijk)-b21(ijk)) - s2*(b20(ijk)+a21(ijk))
     2001                         d(jc+j) = s2*(a20(ijk)-b21(ijk)) + c2*(b20(ijk)+a21(ijk))
     2002                         c(je+j) = c4*(a20(ijk)+b21(ijk)) - s4*(b20(ijk)-a21(ijk))
     2003                         d(je+j) = s4*(a20(ijk)+b21(ijk)) + c4*(b20(ijk)-a21(ijk))
     2004
     2005                         a11(ijk) = (a(ie+i)-a(ib+i)) + (a(ic+i)-a(if+i))
     2006                         b11(ijk) = (b(ie+i)+b(ib+i)) - (b(ic+i)+b(if+i))
     2007                         a20(ijk) = (a(ia+i)-a(id+i)) - 0.5_wp*a11(ijk)
     2008                         a21(ijk) = sin60*((a(ie+i)-a(ib+i))-(a(ic+i)-a(if+i)))
     2009                         b20(ijk) = (b(ia+i)+b(id+i)) + 0.5_wp*b11(ijk)
     2010                         b21(ijk) = sin60*((b(ie+i)+b(ib+i))+(b(ic+i)+b(if+i)))
     2011
     2012                         c(jd+j) = c3*((a(ia+i)-a(id+i))+a11(ijk)) - s3*((b(ia+i)+b(id+i))-b11(ijk))
     2013                         d(jd+j) = s3*((a(ia+i)-a(id+i))+a11(ijk)) + c3*((b(ia+i)+b(id+i))-b11(ijk))
     2014                         c(jb+j) = c1*(a20(ijk)-b21(ijk)) - s1*(b20(ijk)-a21(ijk))
     2015                         d(jb+j) = s1*(a20(ijk)-b21(ijk)) + c1*(b20(ijk)-a21(ijk))
     2016                         c(jf+j) = c5*(a20(ijk)+b21(ijk)) - s5*(b20(ijk)+a21(ijk))
     2017                         d(jf+j) = s5*(a20(ijk)+b21(ijk)) + c5*(b20(ijk)+a21(ijk))
     2018
     2019                         i = i + inc3
     2020                         j = j + inc4
     2021                      ENDDO
     2022                      ibase = ibase + inc1
     2023                      jbase = jbase + inc2
     2024                   ENDDO
     2025                   ia = ia + iink
     2026                   ib = ib + iink
     2027                   ic = ic + iink
     2028                   id = id - iink
     2029                   ie = ie - iink
     2030                   if = if - iink
     2031                   jbase = jbase + jump
     2032                ENDDO
     2033
     2034                IF ( ic > id )  RETURN
     2035
     2036             ENDIF
     2037
     2038             ibase = 0
     2039             DO  l = 1, la
     2040                i = ibase
     2041                j = jbase
     2042                !DIR$ IVDEP
     2043                !CDIR NODEP
     2044                !OCL NOVREC
     2045                DO  ijk = 1, lot
     2046                   c(ja+j) = a(ib+i) + (a(ia+i)+a(ic+i))
     2047                   c(jd+j) = b(ib+i) - (b(ia+i)+b(ic+i))
     2048                   c(jb+j) = (sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i))
     2049                   c(jf+j) = -(sin60*(a(ia+i)-a(ic+i))) - (0.5_wp*(b(ia+i)+b(ic+i))+b(ib+i))
     2050                   c(jc+j) = sin60*(b(ic+i)-b(ia+i)) + (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i))
     2051                   c(je+j) = sin60*(b(ic+i)-b(ia+i)) - (0.5_wp*(a(ia+i)+a(ic+i))-a(ib+i))
     2052                   i = i + inc3
     2053                   j = j + inc4
     2054                ENDDO
     2055                ibase = ibase + inc1
     2056                jbase = jbase + inc2
     2057             ENDDO
     2058
     2059          ELSE
     2060
     2061             ssin60 = 2.0_wp * sin60
     2062             DO  l = 1, la
     2063                i = ibase
     2064                j = jbase
     2065                !DIR$ IVDEP
     2066                !CDIR NODEP
     2067                !OCL NOVREC
     2068                DO  ijk = 1, lot
     2069                   c(ja+j) = (2.0_wp*(a(ia+i)+a(id+i))) + (2.0_wp*(a(ib+i)+a(ic+i)))
     2070                   c(jd+j) = (2.0_wp*(a(ia+i)-a(id+i))) - (2.0_wp*(a(ib+i)-a(ic+i)))
     2071                   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)))
     2072                   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)))
     2073                   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)))
     2074                   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)))
     2075                   i = i + inc3
     2076                   j = j + inc4
     2077                ENDDO
     2078                ibase = ibase + inc1
     2079                jbase = jbase + inc2
     2080             ENDDO
     2081
     2082          ENDIF
     2083
     2084!
     2085!--    Coding for factor 8
     2086       CASE ( 6 )
     2087
     2088          IF ( la /= m )  THEN
     2089             ierr = 3
     2090             RETURN
     2091          ENDIF
     2092
     2093          ia = 1
     2094          ib = ia + la*inc1
     2095          ic = ib + 2*la*inc1
     2096          id = ic + 2*la*inc1
     2097          ie = id + 2*la*inc1
     2098          ja = 1
     2099          jb = ja + jink
     2100          jc = jb + jink
     2101          jd = jc + jink
     2102          je = jd + jink
     2103          jf = je + jink
     2104          jg = jf + jink
     2105          jh = jg + jink
     2106          ssin45 = SQRT( 2.0_wp )
     2107
     2108          DO  l = 1, la
     2109             i = ibase
     2110             j = jbase
     2111             !DIR$ IVDEP
     2112             !CDIR NODEP
     2113             !OCL NOVREC
     2114             DO  ijk = 1, lot
     2115                c(ja+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))+(a(ib+i)+a(id+i)))
     2116                c(je+j) = 2.0_wp*(((a(ia+i)+a(ie+i))+a(ic+i))-(a(ib+i)+a(id+i)))
     2117                c(jc+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))-(b(ib+i)-b(id+i)))
     2118                c(jg+j) = 2.0_wp*(((a(ia+i)+a(ie+i))-a(ic+i))+(b(ib+i)-b(id+i)))
     2119                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)))
     2120                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)))
     2121                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)))
     2122                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)))
     2123                i = i + inc3
     2124                j = j + inc4
     2125             ENDDO
     2126             ibase = ibase + inc1
     2127             jbase = jbase + inc2
     2128          ENDDO
     2129
     2130    END SELECT
     2131
     2132 END SUBROUTINE rpassm
    20812133
    20822134!------------------------------------------------------------------------------!
     
    20882140!> functins required by fft99 & fft991cy
    20892141!------------------------------------------------------------------------------!
    2090   SUBROUTINE set99(trigs,ifax,n)
    2091 
    2092 
    2093     USE control_parameters,                                                    &
     2142 SUBROUTINE set99( trigs, ifax, n )
     2143
     2144
     2145    USE control_parameters,                                                                        &
    20942146        ONLY:  message_string
    20952147
     
    20982150    IMPLICIT NONE
    20992151
    2100     !  Scalar arguments
    21012152    INTEGER(iwp) ::  n !<
    21022153
    2103     !  Array arguments
    21042154    INTEGER(iwp) ::  ifax(*)  !<
    21052155    REAL(wp)     ::  trigs(*) !<
    21062156
    21072157
    2108     !  Local scalars:
    21092158    REAL(wp) ::  angle    !<
    21102159    REAL(wp) ::  del      !<
     
    21192168    INTEGER(iwp) ::  nu   !<
    21202169
    2121     !  Local arrays:
    21222170    INTEGER(iwp) ::  jfax(10) !<
    21232171    INTEGER(iwp) ::  lfax(7)  !<
    21242172
    2125     !  Intrinsic functions
    2126 !    INTRINSIC ASIN, COS, MOD, REAL, SIN
    2127 
    2128     !  Data statements
    2129     DATA lfax/6, 8, 5, 4, 3, 2, 1/
    2130 
    2131 
    2132     !  Executable statements
     2173    DATA  lfax / 6, 8, 5, 4, 3, 2, 1 /
     2174
     2175
    21332176    ixxx = 1
    21342177
    2135     del = 4.0_wp*ASIN(1.0_wp)/REAL(n,KIND=wp)
     2178    del = 4.0_wp * ASIN( 1.0_wp ) / REAL( n, KIND=wp )
    21362179    nil = 0
    21372180    nhl = (n/2) - 1
    2138     DO k = nil, nhl
    2139        angle = REAL(k,KIND=wp)*del
    2140        trigs(2*k+1) = COS(angle)
    2141        trigs(2*k+2) = SIN(angle)
    2142     END DO
    2143 
    2144     ! Find factors of n (8,6,5,4,3,2; only one 8 allowed)
    2145     ! Look for sixes first, store factors in descending order
    2146     nu = n
     2181    DO  k = nil, nhl
     2182       angle = REAL( k, KIND=wp ) * del
     2183       trigs(2*k+1) = COS( angle )
     2184       trigs(2*k+2) = SIN( angle )
     2185    ENDDO
     2186
     2187!
     2188!-- Find factors of n (8,6,5,4,3,2; only one 8 allowed).
     2189!-- Look for sixes first, store factors in descending order.
     2190    nu   = n
    21472191    ifac = 6
    2148     k = 0
    2149     l = 1
    2150 10  CONTINUE
    2151     IF (MOD(nu,ifac)/=0) GO TO 30
    2152     k = k + 1
    2153     jfax(k) = ifac
    2154     IF (ifac/=8) GO TO 20
    2155     IF (k==1) GO TO 20
    2156     jfax(1) = 8
    2157     jfax(k) = 6
    2158 20  CONTINUE
    2159     nu = nu/ifac
    2160     IF (nu==1) GO TO 40
    2161     IF (ifac/=8) GO TO 10
    2162 30  CONTINUE
    2163     l = l + 1
    2164     ifac = lfax(l)
    2165     IF (ifac>1) GO TO 10
    2166 
    2167 !    WRITE (nout,'(A,I4,A)') ' n =',n,' - Contains illegal factors'
    2168     message_string = 'number of gridpoints along x or/and y ' //               &
    2169                      'contain illegal  factors' //                             &
    2170                      '&only factors 2, 3, 5 are allowed'
    2171     CALL message( 'temperton_fft', 'PA0311', 1, 2, 0, 6, 0 )
    2172 
    2173     RETURN
    2174 
    2175     ! Now reverse order of factors
    2176 40  CONTINUE
     2192    k    = 0
     2193    l    = 1
     2194
     2195    DO
     2196
     2197       IF ( MOD( nu, ifac ) == 0 )  THEN
     2198          k = k + 1
     2199          jfax(k) = ifac
     2200          IF ( ifac == 8  .AND.  k /= 1 )  THEN
     2201             jfax(1) = 8
     2202             jfax(k) = 6
     2203          ENDIF
     2204          nu = nu/ifac
     2205          IF ( nu == 1   )  EXIT
     2206          IF ( ifac /= 8 )  CYCLE
     2207       ENDIF
     2208
     2209       l = l + 1
     2210       ifac = lfax(l)
     2211       IF (ifac < 2 )  THEN
     2212          message_string = 'number of gridpoints along x or/and y ' //                             &
     2213                           'contain illegal  factors' //                                           &
     2214                           '&only factors 2, 3, 5 are allowed'
     2215          CALL message( 'temperton_fft', 'PA0311', 1, 2, 0, 6, 0 )
     2216          RETURN
     2217       ENDIF
     2218
     2219    ENDDO
     2220!
     2221!-- Now reverse order of factors
    21772222    nfax = k
    21782223    ifax(1) = nfax
    2179     DO i = 1, nfax
     2224    DO  i = 1, nfax
    21802225       ifax(nfax+2-i) = jfax(i)
    2181     END DO
     2226    ENDDO
    21822227    ifax(10) = n
    2183     RETURN
    2184   END SUBROUTINE set99
     2228
     2229 END SUBROUTINE set99
    21852230
    21862231 END MODULE temperton_fft
  • palm/trunk/SOURCE/wind_turbine_model_mod.f90

    r3685 r3725  
    2626! -----------------
    2727! $Id$
     28! unused variables removed
     29!
     30! 3685 2019-01-21 01:02:11Z knoop
    2831! Some interface calls moved to module_interface + cleanup
    2932!
     
    237240    REAL(wp), DIMENSION(1:100) ::  rnac             = 0.0_wp  !< nacelle diameter [m]
    238241    REAL(wp), DIMENSION(1:100) ::  rr              = 63.0_wp  !< rotor radius [m]
    239     REAL(wp), DIMENSION(1:100) ::  turb_cd_nacelle = 0.85_wp  !< drag coefficient for nacelle
     242!    REAL(wp), DIMENSION(1:100) ::  turb_cd_nacelle = 0.85_wp  !< drag coefficient for nacelle
    240243    REAL(wp), DIMENSION(1:100) ::  turb_cd_tower    = 1.2_wp  !< drag coefficient for tower
    241244
     
    306309    REAL(wp) ::  eps_min          !<
    307310    REAL(wp) ::  eps_min2         !<
    308     REAL(wp) ::  sqrt_arg         !<
    309311
    310312!
     
    519521       IMPLICIT NONE
    520522       
    521        INTEGER(iwp) ::  ierrn       !<
    522 
    523523       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file
    524524
     
    531531                                  rnac, rr, segment_length, segment_width,     &
    532532                                  slope2, speed_control, tilt, time_turbine_on,&
    533                                   turb_cd_nacelle, turb_cd_tower, pitch_rate,  &
     533                                  turb_cd_tower, pitch_rate,                   &
    534534                                  yaw_control, yaw_speed, tl_cor
     535!                                  , turb_cd_nacelle
    535536                                 
    536537       NAMELIST /wind_turbine_parameters/                                      &
     
    543544                                  rnac, rr, segment_length, segment_width,     &
    544545                                  slope2, speed_control, tilt, time_turbine_on,&
    545                                   turb_cd_nacelle, turb_cd_tower, pitch_rate,  &
     546                                  turb_cd_tower, pitch_rate,                   &
    546547                                  yaw_control, yaw_speed, tl_cor
     548!                                  , turb_cd_nacelle
    547549!
    548550!--    Try to find wind turbine model package
     
    980982       INTEGER(iwp) ::  tower_n      !<
    981983       INTEGER(iwp) ::  tower_s      !<
    982 !
    983 !--    Help variables for the calulaction of the nacelle drag
    984        INTEGER(iwp) ::  i_ip         !<
    985        INTEGER(iwp) ::  i_ipg        !<
    986        
    987        REAL(wp) ::  yvalue               
    988        REAL(wp) ::  dy_int           !<
    989        REAL(wp) ::  dz_int           !<
    990984       
    991985       REAL(wp), DIMENSION(:,:), ALLOCATABLE :: circle_points  !<
Note: See TracChangeset for help on using the changeset viewer.