text:
subroutine z_predict_nhfull(j         ,nmmaxj    ,nmmax     ,kmax      , &
                          & mmax      ,nmax      ,nfltyp    ,nst       , &
                          & nocol     ,norow     ,nsrc      ,dismmt    ,irocol    , &
                          & mnksrc    ,kfu       ,kfv       ,kfs       ,            &
                          & kspu      ,kspv      ,kadu      ,kadv      ,kcs       , &
                          & kcu       ,kcv       ,kfsmn0    ,kfsmx0    , &
                          & kfumn0    ,kfumx0    ,kfvmn0    , &
                          & kfvmx0    ,kfuz0     ,kfvz0     ,kfsz0     ,kcu45     , &
                          & kcv45     ,kcscut    ,s0        ,w0        , &
                          & u0        ,u1        ,v0        ,v1        ,hu        , &
                          & hv        ,thick     ,umean     ,ubrlsu    ,ubrlsv    , &
                          & vmean     ,dpu       ,dpv       ,dps       ,dzu0      , &
                          & dzv0      ,dzs0      ,qxk       ,qyk       ,circ2d    , &
                          & circ3d    ,drhodx    ,drhody    ,disch     ,umdis     , &
                          & vmdis     ,wsu       ,wsv       ,dfu       ,dfv       , &
                          & deltau    ,deltav    ,tp        ,rlabda    ,wsbodyu   , &
                          & wsbodyv   ,fxw       ,fyw       ,gud       ,gvd       , &
                          & guu       ,guv       ,gvv       ,gvu       ,guz       , &
                          & gvz       ,gsqs      ,gsqiu     ,gsqiv     ,taubpu    , &
                          & taubpv    ,taubsu    ,taubsv    ,vicuv     ,vnu2d     , &
                          & vicww     ,rxx       ,rxy       ,ryy       ,windu     , &
                          & windv     ,patm      ,fcorio    ,tgfsep    ,wrkb1     , &
                          & wrkb2     ,wrkb3     ,wrkb4     ,wrkb5     ,wrkb6     , &
                          & wrkb7     ,wrkb8     ,wrkb9     ,wrkb10    ,wrkb11    , &
                          & wrkb12    ,wrkb13    ,wrkb14    , &
                          & zk        ,p0        ,crbc      , &
                          & pship     ,diapl     ,rnpl      ,cfurou    ,cfvrou    , &
                          & precip    ,gdp       )

!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation version 3.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program.  If not, see <http://www.gnu.org/licenses/>.
!
!  contact: delft3d.support@deltares.nl
!  Stichting Deltares
!  P.O. Box 177
!  2600 MH Delft, The Netherlands
!
!  All indications and logos of, and references to, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: z_predict_nhfull.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/kernel/src/compute/z_predict_nhfull.f90 $
!!--description-----------------------------------------------------------------
!
!    Function: ADI performs one time step of the Alternating
!              Direction Implicit (ADI) method
! Method used: A.D.I. method is used.
!              Upwind-approach for wet cross section in shallow
!              areas or if the model area contains structures.
!
!!--pseudo code and references--------------------------------------------------
!
! Written by Sebastian Ullmann
!
!!--declarations----------------------------------------------------------------
    use precision
    use flow2d3d_timers
    !
    use globaldata
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    real(fp)     , dimension(:,:)        , pointer :: ustokes
    real(fp)     , dimension(:,:)        , pointer :: vstokes
    include 'flow_steps_f.inc'
!
! Global variables
!
    integer                                               :: j       !!  Begin pointer for arrays which have been transformed into 1D arrays. Due to the shift in the 2nd (M-) index, J = -2*NMAX + 1
    integer                                               :: kmax    !  Description and declaration in esm_alloc_int.f90
    integer                                               :: mmax    !  Description and declaration in esm_alloc_int.f90
    integer                                               :: nfltyp  !  Description and declaration in esm_alloc_int.f90
    integer                                               :: nmax    !  Description and declaration in esm_alloc_int.f90
    integer                                               :: nmmax   !  Description and declaration in dimens.igs
    integer                                               :: nmmaxj  !  Description and declaration in dimens.igs
    integer                                               :: nocol   !  Description and declaration in esm_alloc_int.f90
    integer                                               :: norow   !  Description and declaration in esm_alloc_int.f90
    integer                                               :: nsrc    !  Description and declaration in esm_alloc_int.f90
    integer                                               :: nst     !!  Time step number
    integer, dimension(5, norow + nocol)                  :: irocol  !  Description and declaration in esm_alloc_int.f90
    integer, dimension(7, nsrc)                           :: mnksrc  !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)             :: kcs     !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)             :: kcu     !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)             :: kcv     !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)             :: kfs     !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)             :: kfsmn0  !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)             :: kfsmx0  !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)             :: kfu     !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)             :: kfumn0  !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)             :: kfumx0  !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)             :: kfv     !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)             :: kfvmn0  !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub)             :: kfvmx0  !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax)     :: kspu    !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax)     :: kspv    !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub,   kmax)     :: kadu    !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub,   kmax)     :: kadv    !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub, kmax)       :: kcscut  !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub, kmax)       :: kcu45   !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub, kmax)       :: kcv45   !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub, kmax)       :: kfsz0   !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub, kmax)       :: kfuz0   !  Description and declaration in esm_alloc_int.f90
    integer, dimension(gdp%d%nmlb:gdp%d%nmub, kmax)       :: kfvz0   !  Description and declaration in esm_alloc_int.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, 3)         :: cfurou  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, 3)         :: cfvrou  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(12, norow+nocol)                  :: crbc    !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(4, norow + nocol)                 :: circ2d  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: deltau  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: deltav  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: dfu     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: dfv     !  Description and declaration in esm_alloc_real.f90
    real(prec), dimension(gdp%d%nmlb:gdp%d%nmub)          :: dps     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: dpu     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: dpv     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: fcorio  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: fxw     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: fyw     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: gsqiu   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: gsqiv   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: gsqs    !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: gud     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: guu     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: guv     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: guz     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: gvd     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: gvu     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: gvv     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: gvz     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: hu      !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: hv      !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: patm    !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: precip  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: pship   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: rlabda  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: s0      !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: taubpu  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: taubpv  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: taubsu  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: taubsv  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: tgfsep  !!  Water elev. induced by tide gen.force. Internal work array WRKB17 used
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: tp      !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: umean   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: vmean   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: vnu2d   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: windu   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: windv   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: wsu     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: wsv     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: wsbodyu !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub)            :: wsbodyv !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax)    :: vicww   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax)    :: w0      !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax + 2)  :: vicuv   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: diapl   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: dzs0    !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: dzu0    !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: dzv0    !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: p0      !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: qxk     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: qyk     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: rnpl    !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: rxx     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: rxy     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: ryy     !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: drhodx  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: drhody  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: u0      !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: u1      !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: ubrlsu  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: ubrlsv  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: v0      !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: v1      !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: wrkb1   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: wrkb2   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: wrkb3   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: wrkb4   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: wrkb5   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: wrkb6   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: wrkb7   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: wrkb8   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: wrkb9   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: wrkb10  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: wrkb11  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: wrkb12  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: wrkb13  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(gdp%d%nmlb:gdp%d%nmub, kmax)      :: wrkb14  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(kmax)                             :: thick   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(0:kmax)                           :: zk
    real(fp), dimension(kmax, 2, norow + nocol)           :: circ3d  !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nsrc)                             :: disch   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nsrc)                             :: umdis   !  Description and declaration in esm_alloc_real.f90
    real(fp), dimension(nsrc)                             :: vmdis   !  Description and declaration in esm_alloc_real.f90
    character(1), dimension(nsrc)                         :: dismmt  !  Description and declaration in esm_alloc_char.f90
!
! Local variables
!
    integer :: idummy
    integer :: icx
    integer :: icy
    integer :: idry
    integer :: nhystp
    integer :: nmaxddb
    logical :: flood   ! Flag for activating flooding part of checku subroutine
!
!! executable statements -------------------------------------------------------
!
    ustokes             => gdp%gdtrisol%ustokes
    vstokes             => gdp%gdtrisol%vstokes
    !
    nmaxddb = nmax + 2*gdp%d%ddbound
    !
    ! Calculate HU and set KFU = 0 for HU < HTRSH (.5*DRYFLC)
    ! Calling checku is necessary because the hu calculation
    ! is performed on the full computational domain (kcu=1)
    ! only the wet points (kfu=1) and flooding is allowed.
    !
    idry  = 0
    flood = .true.
    icx   = nmaxddb
    icy   = 1
    call z_checku(j         ,nmmaxj    ,nmmax     ,icx       ,kmax      , &
                & flood     ,kfu       ,kcs       ,kcu       ,kspu      , &
                & kfumn0    ,kfumx0    ,hu        ,s0        ,dpu       , &
                & dps       ,umean     ,kfuz0     ,kfsmn0    ,kfsmx0    , &
                & u0        ,dzu0      ,zk        ,gdp       )
    !
    ! Calculate HV and set KFV = 0 for HV < HTRSH (.5*DRYFLC)
    ! Calling checku is necessary because the hv calculation
    ! is performed on the full computational domain (kcv=1) instead of
    ! only the wet points (kfv=1) and flooding is allowed.
    !
    idry  = 0
    flood = .true.
    icx   = 1
    icy   = nmaxddb
    call z_checku(j         ,nmmaxj    ,nmmax     ,icx       ,kmax      , &
                & flood     ,kfv       ,kcs       ,kcv       ,kspv      , &
                & kfvmn0    ,kfvmx0    ,hv        ,s0        ,dpv       , &
                & dps       ,vmean     ,kfvz0     ,kfsmn0    ,kfsmx0    , &
                & v0        ,dzv0      ,zk        ,gdp       )
    !
    ! Computation of U1, i.e. evaluate momentum equation for one timest
    !     calculate hu and set kfu = 0 for hu < htrsh (.5*dryflc)
    !                          kfu = 1 for hu > trsh
    !
    call timer_start(timer_uzd, gdp)
    gdp%dd%uzditer = 0
    icx            = nmaxddb
    icy            = 1
    call z_uzd(j         ,nmmaxj    ,nmmax     ,kmax      ,icx       , &
             & icy       ,nsrc      ,kcs       ,kcu45     ,kcscut    , &
             & kcu       ,kfu       ,kfuz0     ,kfumn0    ,kfumx0    , &
             & kfv       ,kfvz0     ,kfvmn0    ,kfvmx0    ,dzv0      , &
             & kfs       ,kfsz0     ,kfsmn0    ,kfsmx0    , &
             & u0        ,v0        ,w0        ,hu        ,hv        ,dzu0      ,dzs0      , &
             & guu       ,gvv       ,gvu       ,guv       ,gsqs      , &
             & gud       ,gvd       ,guz       ,gvz       ,gsqiu     , &
             & disch     ,umdis     ,kspu      ,mnksrc    ,dismmt    , &
             & wrkb1     ,wrkb2     ,wrkb3     ,wrkb4     ,wrkb5     , &
             & wrkb6     ,wrkb7     ,wrkb8     ,wrkb9     ,wrkb10    , &
             & wrkb11    ,wrkb12    ,wrkb13    ,wrkb14    ,circ2d(1, 1) , circ3d(1, 1, 1) , &
             & vicuv     ,vnu2d     ,vicww     ,tgfsep    ,dps       , &
             & dfu       ,deltau    ,tp        ,rlabda    ,fxw       ,wsbodyu   , &
             & drhodx    ,wsu       ,taubpu    ,taubsu    ,rxx       , &
             & rxy       ,windu     ,patm      ,fcorio    ,p0        , &
             & ubrlsu    ,pship     ,diapl     ,rnpl      ,cfurou    , &
             & u1        ,s0        ,dpu       ,qxk       ,qyk       , &
             & norow     ,nocol     ,irocol(1, 1)         ,nst       ,umean     , &
             & crbc(1,1) ,ustokes   ,gdp       )
    !
    ! Computation of V1, i.e. evaluate momentum equation for one half timest
    !     calculate hv and set kfv = 0 for hv < htrsh (.5*dryflc)
    !                              = 1 for hv > trsh
    !
    gdp%dd%uzditer = 0
    icx            = 1
    icy            = nmaxddb
    !
    call z_uzd(j         ,nmmaxj    ,nmmax     ,kmax      ,icx       , &
             & icy       ,nsrc      ,kcs       ,kcv45     ,kcscut    , &
             & kcv       ,kfv       ,kfvz0     ,kfvmn0    ,kfvmx0    , &
             & kfu       ,kfuz0     ,kfumn0    ,kfumx0    ,dzu0      , &
             & kfs       ,kfsz0     ,kfsmn0    ,kfsmx0    , &
             & v0        ,u0        ,w0        ,hv        ,hu        ,dzv0      ,dzs0      , &
             & gvv       ,guu       ,guv       ,gvu       ,gsqs      , &
             & gvd       ,gud       ,gvz       ,guz       ,gsqiv     , &
             & disch     ,vmdis     ,kspv      ,mnksrc    ,dismmt    , &
             & wrkb1     ,wrkb2     ,wrkb3     ,wrkb4     ,wrkb5     , &
             & wrkb6     ,wrkb7     ,wrkb8     ,wrkb9     ,wrkb10    , &
             & wrkb11    ,wrkb12    ,wrkb13    ,wrkb14    ,circ2d(1, norow+1) , circ3d(1, 1, norow+1) , &
             & vicuv     ,vnu2d     ,vicww     ,tgfsep    ,dps       , &
             & dfv       ,deltav    ,tp        ,rlabda    ,fyw       ,wsbodyv   , &
             & drhody    ,wsv       ,taubpv    ,taubsv    ,ryy       , &
             & rxy       ,windv     ,patm      ,fcorio    ,p0        , &
             & ubrlsv    ,pship     ,diapl     ,rnpl      ,cfvrou    , &
             & v1        ,s0        ,dpv       ,qyk       ,qxk       , &
             & nocol     ,norow     ,irocol(1, norow + 1) ,nst       ,vmean     , &
             & crbc(1,norow + 1)    ,vstokes   ,gdp       )
    call timer_stop(timer_uzd, gdp)
    !
    ! DD code added:
    !
    !
    ! Synchronize on Dry Point
    if (nfltyp == 0) then
       nhystp = nxtdry(d3dflow_check_adi_dry, 0, gdp)
    else
       nhystp = nxtdry(d3dflow_check_adi_dry, idry, gdp)
    endif
end subroutine z_predict_nhfull

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


! See https://gcc.gnu.org/onlinedocs/gfortran/Interoperability-with-C.html
module test_bind_c
  use iso_c_binding

  ! test derived type
  type, bind(c) :: atype
    real(c_double) :: a
  end type atype
  ! this is
  ! struct {
  !   double a;
  ! } atype;
  ! in C language

  ! test interoperable global variable
  integer(C_INT), bind(C, name="_MyProject_flags") :: global_flag

  ! test Interoperable Subroutines and Functions
  interface
    subroutine print_c(string) bind(C, name="print_C")
      use iso_c_binding, only: c_char
      character(kind=c_char) :: string(*)
    end subroutine print_c
  end interface
end module test_bind_c


text:
!	Standard error for entropy per item under complex sampling.
!	This version is designed for random sampling of psu's within strata.
!	Sampling is with replacement.
!	npsu is the number of psus's per stratum.
!	nstratum is the number of strata.
!	stratum is observation stratum.
!	psu is observation psu.
!	stratify is .true. for stratified sampling.
!	usepsu is .true. for primary sampling units within strata.

!	obsweight is  observation weights.
!	prob is array of observed probabilities.
!	totalitems is weighted sum of items presented.
real(kind=8) function sdent_complex(npsu,nstratum,psu,stratum,stratify,usepsu,obsweight,prob,totalitems)
implicit none
interface

!	The following description is for the function's basic purpose.
!	Its use in this subroutine is a bit different.
!	Estimate the information matrix for a complex sample.
!	This version is designed for random sampling of psu's within strata.
!	Sampling is with replacement.
!	npsu is the number of psus's per stratum.
!	nstratum is the number of strata.
!	psu indicates the psu of an observation.
!	stratum indicates the stratum of an observation.
!	stratify is .true. for stratified sampling.
!	usepsu is .true. for primary sampling units within strata.
!	grads is the array of observation gradients.
!	obsweight contains observation weights.

    function information_complex(npsu,nstratum,psu,stratum,stratify,usepsu,grads,obsweight)
        implicit none
        integer,intent(in)::nstratum
        integer,intent(in),optional::npsu(:),psu(:),stratum(:)
        logical,intent(in)::stratify,usepsu
        real(kind=8),intent(in)::grads(:,:),obsweight(:)
        real(kind=8)::information_complex(size(grads,1),size(grads,1))
    end function information_complex
end interface

integer,intent(in)::nstratum
integer,intent(in),optional::npsu(:),psu(:),stratum(:)
logical,intent(in)::stratify,usepsu
real(kind=8),intent(in)::obsweight(:),prob(:),totalitems

!	al indicates allocation error.
!   nobs is the number of observations.
!	obs is an observation counter.

integer::al,nobs,obs
!   Variance goes to covobs.
!	Log probabilities.
real(kind=8)::covobs(1,1)
real(kind=8),allocatable::obsmat(:,:)
nobs=size(prob)
allocate(obsmat(1,nobs),stat=al)
if(al/=0) stop "Unable to allocate arrays for asymptotic variances computations for complex sampling."
obsmat=0.0_8
do obs=1,nobs
    if(prob(obs)>0.0_8)obsmat(1,obs)=log(prob(obs))
end do
covobs=information_complex(npsu,nstratum,psu,stratum,stratify,usepsu,&
    obsmat(1:1,:),obsweight)
sdent_complex=sqrt(covobs(1,1))/totalitems
deallocate(obsmat)

return
end function sdent_complex

max_stars_repo_name: EducationalTestingService/MIRT

max_issues_repo_name: EducationalTestingService/MIRT

max_forks_repo_name: EducationalTestingService/MIRT


!-------------------------------------------------------------------------------
! Copyright (c) 2016 The University of Tokyo
! This software is released under the MIT License, see LICENSE.txt
!-------------------------------------------------------------------------------
!> \brief Coupling Interface

module hecmw_couple_info_f

  use hecmw_util

  implicit none
  private
  public :: hecmw_couple_get_unit_id
  public :: hecmw_couple_is_member
  public :: hecmw_couple_is_unit_member
  public :: hecmw_couple_is_unit_member_u
  public :: hecmw_couple_is_root
  public :: hecmw_couple_is_unit_root
  public :: hecmw_couple_is_unit_root_u
  public :: hecmw_intercomm_get_size
  public :: hecmw_intracomm_get_size
  public :: hecmw_intracomm_get_size_u
  public :: hecmw_intercomm_get_rank
  public :: hecmw_intracomm_get_rank
  public :: hecmw_intracomm_get_rank_u
  public :: hecmw_intercomm_get_comm
  public :: hecmw_intracomm_get_comm
  public :: hecmw_intracomm_get_comm_u
  public :: hecmw_intercomm_get_group
  public :: hecmw_intracomm_get_group
  public :: hecmw_intracomm_get_group_u

contains

!...
!...
!...
subroutine hecmw_couple_get_unit_id(boundary_id, unit_specifier, unit_id)

  character(len=HECMW_NAME_LEN), intent(in)  :: boundary_id
  integer(kind=kint),            intent(in)  :: unit_specifier
  character(len=HECMW_NAME_LEN), intent(out) :: unit_id
  integer(kind=kint)                         :: ierr

  call hecmw_couple_get_unit_id_if(boundary_id, unit_specifier, unit_id, ierr)
  if(ierr /= 0) call hecmw_abort(hecmw_comm_get_comm())

end subroutine hecmw_couple_get_unit_id

!...
!...
!...
function hecmw_couple_is_member(boundary_id) result (is_member)

  character(len=HECMW_NAME_LEN), intent(in) :: boundary_id
  integer(kind=kint)                        :: is_member

  call hecmw_couple_is_memb_if(boundary_id, is_member)

end function hecmw_couple_is_member

!...
!...
!...
function hecmw_couple_is_unit_member(boundary_id, unit_specifier) result(is_member)

  character(len=HECMW_NAME_LEN), intent(in) :: boundary_id
  integer(kind=kint),            intent(in) :: unit_specifier
  integer(kind=kint)                        :: is_member

  call hecmw_couple_is_unit_memb_if(boundary_id, unit_specifier, is_member)

end function hecmw_couple_is_unit_member

!...
!...
!...
function hecmw_couple_is_unit_member_u(unit_id) result(is_member)

  character(len=HECMW_NAME_LEN), intent(in) :: unit_id
  integer(kind=kint)                        :: is_member

  call hecmw_couple_is_unit_memb_u_if(unit_id, is_member)

end function hecmw_couple_is_unit_member_u

!...
!...
!...
function hecmw_couple_is_root(boundary_id) result(is_root)

  character(len=HECMW_NAME_LEN), intent(in) :: boundary_id
  integer(kind=kint)                        :: is_root

  call hecmw_couple_is_root_if(boundary_id, is_root)

end function hecmw_couple_is_root

!...
!...
!...
function hecmw_couple_is_unit_root(boundary_id, unit_specifier) result(is_root)

  character(len=HECMW_NAME_LEN), intent(in) :: boundary_id
  integer(kind=kint),            intent(in) :: unit_specifier
  integer(kind=kint)                        :: is_root

  call hecmw_couple_is_unit_root_if(boundary_id, unit_specifier, is_root)

end function hecmw_couple_is_unit_root

!...
!...
!...
function hecmw_couple_is_unit_root_u(unit_id) result(is_root)

  character(len=HECMW_NAME_LEN), intent(in) :: unit_id
  integer(kind=kint)                        :: is_root

  call hecmw_couple_is_unit_root_u_if(unit_id, is_root)

end function hecmw_couple_is_unit_root_u

!...
!...
!...
function hecmw_intercomm_get_size(boundary_id) result(psize)

  character(len=HECMW_NAME_LEN), intent(in) :: boundary_id
  integer(kind=kint)                        :: psize

  call hecmw_intercomm_get_size_if(boundary_id, psize)

end function hecmw_intercomm_get_size

!...
!...
!...
function hecmw_intracomm_get_size(boundary_id, unit_specifier) result(psize)

  character(len=HECMW_NAME_LEN), intent(in) :: boundary_id
  integer(kind=kint),            intent(in) :: unit_specifier
  integer(kind=kint)                        :: psize

  call hecmw_intracomm_get_size_if(boundary_id, unit_specifier, psize)

end function hecmw_intracomm_get_size

!...
!...
!...
function hecmw_intracomm_get_size_u(unit_id) result(psize)

  character(len=HECMW_NAME_LEN), intent(in) :: unit_id
  integer(kind=kint)                        :: psize

  call hecmw_intracomm_get_size_u_if(unit_id, psize)

end function hecmw_intracomm_get_size_u

!...
!...
!...
function hecmw_intercomm_get_rank(boundary_id) result(rank)

  character(len=HECMW_NAME_LEN), intent(in) :: boundary_id
  integer(kind=kint)                        :: rank

  call hecmw_intercomm_get_rank_if(boundary_id, rank)

end function hecmw_intercomm_get_rank

!...
!...
!...
function hecmw_intracomm_get_rank(boundary_id, unit_specifier) result(rank)

  character(len=HECMW_NAME_LEN), intent(in) :: boundary_id
  integer(kind=kint),            intent(in) :: unit_specifier
  integer(kind=kint)                        :: rank

  call hecmw_intracomm_get_rank_if(boundary_id, unit_specifier, rank)

end function hecmw_intracomm_get_rank

!...
!...
!...
function hecmw_intracomm_get_rank_u(unit_id) result(rank)

  character(len=HECMW_NAME_LEN), intent(in) :: unit_id
  integer(kind=kint)                        :: rank

  call hecmw_intracomm_get_rank_u_if(unit_id, rank)

end function hecmw_intracomm_get_rank_u

!...
!...
!...
function hecmw_intercomm_get_comm(boundary_id) result(comm)

  character(len=HECMW_NAME_LEN), intent(in) :: boundary_id
  integer(kind=kint)                        :: comm

  call hecmw_intercomm_get_comm_if(boundary_id, comm)

end function hecmw_intercomm_get_comm

!...
!...
!...
function hecmw_intracomm_get_comm(boundary_id, unit_specifier) result(comm)

  character(len=HECMW_NAME_LEN), intent(in) :: boundary_id
  integer(kind=kint),            intent(in) :: unit_specifier
  integer(kind=kint)                        :: comm

  call hecmw_intracomm_get_comm_if(boundary_id, unit_specifier, comm)

end function hecmw_intracomm_get_comm

!...
!...
!...
function hecmw_intracomm_get_comm_u(unit_id) result(comm)

  character(len=HECMW_NAME_LEN), intent(in) :: unit_id
  integer(kind=kint)                        :: comm

  call hecmw_intracomm_get_comm_u_if(unit_id, comm)

end function hecmw_intracomm_get_comm_u

!...
!...
!...
function hecmw_intercomm_get_group(boundary_id) result(group)

  character(len=HECMW_NAME_LEN), intent(in) :: boundary_id
  integer(kind=kint)                        :: group

  call hecmw_intercomm_get_group_if(boundary_id, group)

end function hecmw_intercomm_get_group

!...
!...
!...
function hecmw_intracomm_get_group(boundary_id, unit_specifier) result(group)

  character(len=HECMW_NAME_LEN), intent(in) :: boundary_id
  integer(kind=kint),            intent(in) :: unit_specifier
  integer(kind=kint)                        :: group

  call hecmw_intracomm_get_group_if(boundary_id, unit_specifier, group)

end function hecmw_intracomm_get_group

!...
!...
!...
function hecmw_intracomm_get_group_u(boundary_id) result(group)

  character(len=HECMW_NAME_LEN), intent(in) :: boundary_id
  integer(kind=kint)                        :: group

  call hecmw_intracomm_get_group_u_if(boundary_id, group)

end function hecmw_intracomm_get_group_u

end module hecmw_couple_info_f


!*==zpocon.f90  processed by SPAG 7.51RB at 20:08 on  3 Mar 2022
!> \brief \b ZPOCON
!
!  =========== DOCUMENTATION ===========
!
! Online html documentation available at
!            http://www.netlib.org/lapack/explore-html/
!
!> \htmlonly
!> Download ZPOCON + dependencies
!>
!> [TGZ]
!>
!> [ZIP]
!>
!> [TXT]
!> \endhtmlonly
!
!  Definition:
!  ===========
!
!       SUBROUTINE ZPOCON( UPLO, N, A, LDA, ANORM, RCOND, WORK, RWORK,
!                          INFO )
!
!       .. Scalar Arguments ..
!       CHARACTER          UPLO
!       INTEGER            INFO, LDA, N
!       DOUBLE PRECISION   ANORM, RCOND
!       ..
!       .. Array Arguments ..
!       DOUBLE PRECISION   RWORK( * )
!       COMPLEX*16         A( LDA, * ), WORK( * )
!       ..
!
!
!> \par Purpose:
!  =============
!>
!> \verbatim
!>
!> ZPOCON estimates the reciprocal of the condition number (in the
!> 1-norm) of a complex Hermitian positive definite matrix using the
!> Cholesky factorization A = U**H*U or A = L*L**H computed by ZPOTRF.
!>
!> An estimate is obtained for norm(inv(A)), and the reciprocal of the
!> condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
!> \endverbatim
!
!  Arguments:
!  ==========
!
!> \param[in] UPLO
!> \verbatim
!>          UPLO is CHARACTER*1
!>          = 'U':  Upper triangle of A is stored;
!>          = 'L':  Lower triangle of A is stored.
!> \endverbatim
!>
!> \param[in] N
!> \verbatim
!>          N is INTEGER
!>          The order of the matrix A.  N >= 0.
!> \endverbatim
!>
!> \param[in] A
!> \verbatim
!>          A is COMPLEX*16 array, dimension (LDA,N)
!>          The triangular factor U or L from the Cholesky factorization
!>          A = U**H*U or A = L*L**H, as computed by ZPOTRF.
!> \endverbatim
!>
!> \param[in] LDA
!> \verbatim
!>          LDA is INTEGER
!>          The leading dimension of the array A.  LDA >= max(1,N).
!> \endverbatim
!>
!> \param[in] ANORM
!> \verbatim
!>          ANORM is DOUBLE PRECISION
!>          The 1-norm (or infinity-norm) of the Hermitian matrix A.
!> \endverbatim
!>
!> \param[out] RCOND
!> \verbatim
!>          RCOND is DOUBLE PRECISION
!>          The reciprocal of the condition number of the matrix A,
!>          computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
!>          estimate of the 1-norm of inv(A) computed in this routine.
!> \endverbatim
!>
!> \param[out] WORK
!> \verbatim
!>          WORK is COMPLEX*16 array, dimension (2*N)
!> \endverbatim
!>
!> \param[out] RWORK
!> \verbatim
!>          RWORK is DOUBLE PRECISION array, dimension (N)
!> \endverbatim
!>
!> \param[out] INFO
!> \verbatim
!>          INFO is INTEGER
!>          = 0:  successful exit
!>          < 0:  if INFO = -i, the i-th argument had an illegal value
!> \endverbatim
!
!  Authors:
!  ========
!
!> \author Univ. of Tennessee
!> \author Univ. of California Berkeley
!> \author Univ. of Colorado Denver
!> \author NAG Ltd.
!
!> \date December 2016
!
!> \ingroup complex16POcomputational
!
!  =====================================================================
      SUBROUTINE ZPOCON(Uplo,N,A,Lda,Anorm,Rcond,Work,Rwork,Info)
      IMPLICIT NONE
!*--ZPOCON124
!
!  -- LAPACK computational routine (version 3.7.0) --
!  -- LAPACK is a software package provided by Univ. of Tennessee,    --
!  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
!     December 2016
!
!     .. Scalar Arguments ..
      CHARACTER Uplo
      INTEGER Info , Lda , N
      DOUBLE PRECISION Anorm , Rcond
!     ..
!     .. Array Arguments ..
      DOUBLE PRECISION Rwork(*)
      COMPLEX*16 A(Lda,*) , Work(*)
!     ..
!
!  =====================================================================
!
!     .. Parameters ..
      DOUBLE PRECISION ONE , ZERO
      PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
!     ..
!     .. Local Scalars ..
      LOGICAL upper
      CHARACTER normin
      INTEGER ix , kase
      DOUBLE PRECISION ainvnm , scale , scalel , scaleu , smlnum
      COMPLEX*16 zdum
!     ..
!     .. Local Arrays ..
      INTEGER isave(3)
!     ..
!     .. External Functions ..
      LOGICAL LSAME
      INTEGER IZAMAX
      DOUBLE PRECISION DLAMCH
      EXTERNAL LSAME , IZAMAX , DLAMCH
!     ..
!     .. External Subroutines ..
      EXTERNAL XERBLA , ZDRSCL , ZLACN2 , ZLATRS
!     ..
!     .. Intrinsic Functions ..
      INTRINSIC ABS , DBLE , DIMAG , MAX
!     ..
!     .. Statement Functions ..
      DOUBLE PRECISION CABS1
!     ..
!     .. Statement Function definitions ..
      CABS1(zdum) = ABS(DBLE(zdum)) + ABS(DIMAG(zdum))
!     ..
!     .. Executable Statements ..
!
!     Test the input parameters.
!
      Info = 0
      upper = LSAME(Uplo,'U')
      IF ( .NOT.upper .AND. .NOT.LSAME(Uplo,'L') ) THEN
         Info = -1
      ELSEIF ( N<0 ) THEN
         Info = -2
      ELSEIF ( Lda<MAX(1,N) ) THEN
         Info = -4
      ELSEIF ( Anorm<ZERO ) THEN
         Info = -5
      ENDIF
      IF ( Info/=0 ) THEN
         CALL XERBLA('ZPOCON',-Info)
         RETURN
      ENDIF
!
!     Quick return if possible
!
      Rcond = ZERO
      IF ( N==0 ) THEN
         Rcond = ONE
         RETURN
      ELSEIF ( Anorm==ZERO ) THEN
         RETURN
      ENDIF
!
      smlnum = DLAMCH('Safe minimum')
!
!     Estimate the 1-norm of inv(A).
!
      kase = 0
      normin = 'N'
      DO
         CALL ZLACN2(N,Work(N+1),Work,ainvnm,kase,isave)
         IF ( kase/=0 ) THEN
            IF ( upper ) THEN
!
!           Multiply by inv(U**H).
!
               CALL ZLATRS('Upper','Conjugate transpose','Non-unit',    &
     &                     normin,N,A,Lda,Work,scalel,Rwork,Info)
               normin = 'Y'
!
!           Multiply by inv(U).
!
               CALL ZLATRS('Upper','No transpose','Non-unit',normin,N,A,&
     &                     Lda,Work,scaleu,Rwork,Info)
            ELSE
!
!           Multiply by inv(L).
!
               CALL ZLATRS('Lower','No transpose','Non-unit',normin,N,A,&
     &                     Lda,Work,scalel,Rwork,Info)
               normin = 'Y'
!
!           Multiply by inv(L**H).
!
               CALL ZLATRS('Lower','Conjugate transpose','Non-unit',    &
     &                     normin,N,A,Lda,Work,scaleu,Rwork,Info)
            ENDIF
!
!        Multiply by 1/SCALE if doing so will not cause overflow.
!
            scale = scalel*scaleu
            IF ( scale/=ONE ) THEN
               ix = IZAMAX(N,Work,1)
               IF ( scale<CABS1(Work(ix))*smlnum .OR. scale==ZERO ) EXIT
               CALL ZDRSCL(N,scale,Work,1)
            ENDIF
            CYCLE
         ENDIF
!
!     Compute the estimate of the reciprocal condition number.
!
         IF ( ainvnm/=ZERO ) Rcond = (ONE/ainvnm)/Anorm
         EXIT
      ENDDO
!
!
!     End of ZPOCON
!
      END SUBROUTINE ZPOCON


text:
      program dnff

      character*10 cdate
      character*4 cvers
      parameter (cvers = '5.1', cdate = '02/13/1998')

      include '../../D/iounits.inc'
      include 'params1.inc'
      include '../../D/decimate.inc'
      include '../../D/input.inc'
      include 'dnff_seed.inc'

      parameter (nt=nwmx/2,ntot= nwmx*nchmx,nchmx1=nchmx+1,
     &  nmax = 3*nwmx,nwsvmx = nwmx*4+15,nfmax=nwmx/2,nwmx1=nwmx+10
     &       , nchmx2 = 2*nchmx,nchmx21 = nchmx*2+1)

      real x(nchmx,nwmx1),samprate(ndmx),w(nwmx,ndmx)
     &    ,sc(nchmx),wsave(nwsvmx,ndmx),stcor(2),orient(2,nchmx),
     &    xx(ndmx,nchmx1,0:nxxmx),pspecl1(nchmx,nwmx,ndmx),
     &    areg(nchmx,10),rr(nwmx),ar(nwmx),bwo(nwmx),br(nwmx,1)
      real	dt

      complex wrk(nwmx),rnrmt(nchmx,ndmx,nt),cwrk(nchmx,nsmax)
     &   ,cout(nchmx,nsmax)

      integer ix1(nchmx1,nixmx),
     &   next(ndmx),ifirstd(ndmx),nacc(ndmx),
     &    idp(0:ndmx,2,nbadmx),iwrk(nchmx,nsmax),
     &    ids(0:ndmx,2,nbadmx),ndp(0:ndmx),ifirst(ndmx),nspec(ndmx),
     &    nds(0:ndmx),iftype(nfilmax,nchmx),
     &    nfil(nchmx),iset(4,nsmax),ifdir(3,nsmax),
     &    idl(ndmx),ioffc,ixs(nsmax),doy_start,hr_start,min_start,
     &     sec_start,iset_offset(ndmx),time_int,doy_0

        logical lstart(ndmx),lstartd(ndmx),lgood,lg,lclkd
     &  ,lfd(nchmx,ndmx),lwinst(ndmx),lseed

      character*1 ans
      character*50 arg
      character*6 chid(nchmx)
      character*80 afparam(nfilmax,nchmx),cfin(nchmx),cf_station

C   SCRATCH FILE 02.01.98
        real x_scr(nsmax*nfmax,nchmx,2)
        integer iuse_scr(nsmax*nfmax)

       print*,'********************************************************'
       print*,'*                                                      *'
       print*,'*     DNFF_SEED     Version: ',cvers,'   Date: ',cdate,
     &'        *'
       print*,'*                                                      *'
       print*,'********************************************************'
       print*,' '
       bytes = NBYTES
      ioff = 0
      lseed = .false.

c     parse command line options
       nargs = iargc()
       do k = 1, nargs
          call getarg(k,arg)
          if (arg(1:2) .eq. '-b') then
             read (arg(3:3),*) bytes
          elseif(arg(1:2).eq.'-s') then
ccc          Read directly from mini-seed data files.  In this case only one
ccc          run can be done at a time.  To run with seed files:
ccc                    dnff -scf_station
ccc          where "cf_station" is the name of a file containing # of
ccc          channels, and file names for each channel
             lseed = .true.
             cf_station = arg(3:50)
             if(iclong(cf_station,78).eq.0) cf_station = 'cf_station'
          endif

      enddo
      if( lseed) then
         write(6,'(/,''    Using binary mini-SEED files'',/)')
      else
      write(6,'(/,''    Using integer*'',i1,'' as binary input '',/)')
     &      bytes
      endif

c       return here to start processing another data file
1000  continue

c*********************************************************************

c    compute l1 power spectrum for pre-scaling output
c          blank array pspecl1
      do i = 1,ndmx
         nspec(i) = 0
         do k = 1,nchmx
            do j = 1,nwmx
               pspecl1(k,j,i) = 0.0
            enddo
         enddo
      enddo

c*****  get station name from standard input, input files paths from
c       local file paths; open input file, make sp br etc path/file names
      if(lseed) then
         call init_files(cf_station,cfin,nch,doy_start,hr_start,
     &                   min_start,sec_start)
         call init_seed(cfin,nch,dr,sampfreq,doy_0)
         time_int = - sec_start+
     &       60*(-min_start+60*(-hr_start + 24*(doy_0-doy_start)))
         write(0,*) 'Survey Start'
         write(0,*) doy_start,hr_start,min_start,sec_start
         write(0,*) 'doy_0',doy_0
         write(0,*) 'Time Interval',time_int

         msval = 2147483647
         rmsval =  float(msval) - 1000.
      else
         call cininit(msval,rmsval,nch)
      endif
      nchp1 = nch+1

      write(0,*) 'Opening output file',cfout
C     ME  REMOVE SCRATCH FILE
      call outinit(nch,lpack)

c  read system parameter file:  number of channels,
c  electrode line lengths, filter parameters, conversion
c  factors, etc. from file spsta###x
      call getsp(nch0,sampr,sc,nfil,iftype,afparam,
     &   decl,stcor,orient,chid,cda,cdb,lclkd)  !cfsp,

      if(nch0.ne.nch) then
         print*,'error: nch in sp file does not agree with nch in',
     &     ' data file'
         stop
      endif
      if(cfsp(4:11).eq.'standard') then
         print*,'enter declination, station coordinates'
         read(5,*) decl,stcor(1),stcor(2)
      endif

c  set up decimation; decset reads parameters which control decimation
c  from a file; see documentation in decset for further info;
      call decset(nwmx,lfd,idoff,ierr)
      if(ierr.eq.1) then
         print*,'parameters read from cf_decset are inconsistent',
     &   ' with parameters declared in main program'
           stop
      endif

c  compute sampling rate for each decimation level
        nfusemx = 0
        do 1 i = 1,nd
        idl(i) = i + idoff
        samprate(i) = sampr
        nfusemx = max(nfusemx,nfuse(i))
           do 1 j = 1,i
1          samprate(i) = samprate(i)*idec(j)
c  set up decimation filter corrections, conversion of counts
c  to physical units (this routine makes a table of transfer
c  functions to correct for all of these)
      call fcorsu(samprate,nfil,afparam,
     1  iftype,sc,rnrmt,cda,lfd)
c    read in bad records file (gives record ranges which should not
c    be processed (idp) or stacked (ids); see routine for input format
c    file also gives offset to adjust record numbers with
      call bdrcsu(nd,nbadmx,ndp,idp,nds,ids,ioffc)  !cfbr,
ccccEGBERT
ccc    Change call to mk_offset, and replace mk_offst.f with new version
ccc     declare dt (output by mk_offset, and used by pshftsd) as real
           call mk_offset(time_int,iset_offset,sampfreq,dt)
ccccEGBERT
c      write(6,*) 'TIME_INTERVAL',time_int
c      write(6,*) 'ISET_OFFSET',iset_offset

c     set up fft tables
      do id = 1,nd
         call cffti(nwin(id),wsave(1,id))
      enddo

c     initialize various arrays
      do i = 1,nd
         lwinst(i) = .true.
      enddo
      nuse = 0
      nmsmx = nwmx

c    get first block of data and set up pointers at start of processing
c    also, if too large a data gap occurs return to stmt. 90
c    and reset pointers

      lfirst = .true.
85    continue
      if(lseed) then
         call align_SEED(npts,nch,ix1,ierr)
         call rdblk_seed(npts,nch,ix1,ierr)
         if(npts .le. 0 ) then
            write(0,*) 'Terminating after initial read'
            write(0,*) 'npts,ierr',npts,ierr
            stop
         endif
      else
         call rdblk(npts,ix1,ierr,nch)
      endif
      if(ierr.eq.-1) go to 300
      irec1 = ix1(1,1)
         write(6,*) 'Time interval since t0',time_int
         write(6,*) 'First sample number ',irec1

90    call pterst(irec1,ifirst,ifirstd,next,lstart,lstartd)
      write(0,*) 'Setup Complete'

c>>>>>>>>>>>>>>>main loop starts here >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
        ifrec = 0
100     continue

c     call decimation routine
        call dcimte(ix1,npts,xx,nchp1,nxxmx,ifirst,ifirstd,next,lstart
     &    ,lstartd,nacc,rmsval)

      do 200 id = 1,nd
      if(nfuse(id) .gt. 0) then
c     compute number of sets accumulated at decimation level id
c    (will ususally be 1 or 0)
      nsets = (nacc(id) - olap(id) - 2*npwmx(id) - missmx(1,id)/2)
     &                           / ( nwin(id) - olap(id) )
         do 190 i = 1,nsets
c          make set i at decimation level id; return in array x
           call mkset(xx,nchp1,nxxmx,ifirst,x,nwmx,ist,id,irec,lrec,
     &                rmsval,lgood,roff)

            ist = ist + iset_offset(id)
c   need to decide if set contains data segmenta which should not be
c   processed (declared in bad record file brsta###x)
           call badrec(nd,ndp,idp,id,irec,lrec,lg)
           if (lgood.and.lg) then
c       set is "good" (i.e. not too many missing data points; not
c       among the bad record ranges); proceed with processing

c   first difference,demean, window, and fourier transform data
             n = nwin(id)
             nmen=n+npwmx(id)-1
             call demean(x,nch,nmen)
             call frstdif(x,nch,n,lfd,id,areg,npw(1,id),rr,ar,br,bwo)
             call demean(x,nch,n)
             call ps1win(x,nch,n,lwinst(id),w(1,id))
             call fftnp2(n,nch,x,wrk,wsave(1,id))

c   correct for filters (decimation, first differencing) and normalize
c   fourier coefficients to have units of nt/(sqrt(hz))
c   [ uses table of transfer function coefficients stored in complex
c     array rnrmt ]
            nfmid = nfuse(id)
            call filtcor(x,nch,nfmid,id,rnrmt,nd,areg,npw(1,id),n,lfd)
            nfreq = nfmid

CCCCC   WHAT IS CORRECT HERE ??????
ccccEGBERT
ccc ADD THIS CALL ... and include new subroutine pshftsd.f
cc      correct for offsets in actual sampling times

ccc      for seed files need to correct for offsets in actual sampling
ccc       times relative to nominal sampling times
             if(lseed) then
ccc             use seed version
                call pshftsd(x,nfreq,nch,samprate(id),nwin(id),irec)
             else
               call phs_shft(dt,x,nfreq,nch,samprate(id),nwin(id))
             endif

ccccEGBERT

c        if necessary correct for clock drifts
              if(lclkd) then
                 tset = ((irec+lrec)/2. - iclk0)*samprate(1)
                 call cldrft(x,nfreq,nch,samprate(id)
     &              ,nwin(id),cdb,tset)
              endif

c     if necessary correct for offset of begining of set (caused by
c        non-integer set overlap)
              if(abs(roff) .gt. 1.0e-5) then
                 roff = - roff * samprate(1)
                 call cldrft(x,nfreq,nch,samprate(id)
     &              ,nwin(id),1.0,roff)
              endif

c       check to see if set contains segments which should not be stacked
c       if so indicate by making set number negative
              call badrec(nd,nds,ids,id,irec,lrec,lg)
              if( .not. lg) ist = - abs(ist)

c         output fourier coefficients to scratch file f_sta###x
              call freqout(nch,nfreq,x,ifrec,nsmax,nfusemx,
     &            x_scr,iuse_scr)

c         add abs(FCs) to array to compute scalling factor for output
              call l1spec(nch,x,nfreq,pspecl1(1,1,id),nchmx)
              nspec(id) = nspec(id) + 1

c         update array ifdir (which will be used to sort out fourier coefficients
c         at end of program)
              nuse = nuse + 1
              ifdir(1,nuse) = id + idoff
              ifdir(2,nuse) = nfreq
              ifdir(3,nuse) = ist
c     check for overflow of nsmax (maximum # of sets
              if (nuse.eq.nsmax) then
                 write(*,232) nsmax
 232   format(/,'!!!  Number of processed sets reached',/,
     &          '     maximum number of allowed sets !!!',/,
     &          '     The program will finish processing now',/,
     &          '     To process all data change NSMAX =',i6,/,
     &          '     in file PARAMS1.INC !!!')
                 goto 300

               endif
c              print*,'good',irec,lrec,id,ist
c           else
c              print*,'bad',irec,lrec,id,ist
           endif
190     continue
        endif
200     continue
c    get more data and go back to top of main loop
      if(lseed) then
         call rdblk_seed(npts,nch,ix1,ierr)
      else
         call rdblk(npts,ix1,ierr,nch)
      endif

c     ierr = 0 for normal read
        if(ierr.eq.0) then
           lfirst = .false.
           go to 100
        endif

c      ierr = -3 when there is a gap in the record numbers that is so
c large that it is most reasonable to start the decimation filtering over
        if(ierr.eq. -3) then
           if(lseed) then
              lfirst = .true.
              go to 85
           else
              go to 90
           endif
        endif

c     ierr = -2 means the record numbers jump backwards; the program doesn't
c    know what to do about this and stops
ccc     won't happen with seed files
        if(ierr.eq.-2) then
            print*,'error: excess records/ incorrect time mark'
            print*,'stopping here'
        endif

300     continue

c<<<<<<<<<<<<<<<<<<< end of main loop <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      print*,'number of sets = ',nuse

c***********************
c       reorder f_ file by frequency
c*************************

c       close connection to data file
      close(in_unit)
c       reorder and output file of fourier coefficients
c      l1 power spectrum
      do 490 i = 1,nd
      do 490 j = 1,nch
      do 490 k = 1,nfuse(i)
         if(lpack) then
            pspecl1(j,k,i) = pspecl1(j,k,i)/nspec(i)
         else
            pspecl1(j,k,i) = 1000.
         endif
490      continue

ccc    eliminate decimation levels with nfuse(id) = 0
      iid = 0
      do id = 1,nd
         if(nfuse(id) .gt. 0) then
            iid = iid + 1
            idl(iid) = idl(id)
            nwin(iid) = nwin(id)
            samprate(iid) = samprate(id)
         endif
      enddo

c      header for FC file
      irlo = 4*(nch+1)
      call wfhead(nch,iid,nfmax,nwin,samprate,idl
     1   ,chid,orient,decl,stcor,irlo)

c       read in scratch file, reorder pack if requested, output
C     ME 	REMOVE SCRATCH FILE
c      call mkseq(ioscr,iouout,ifdir,iset,nuse,nch,nd,
c     1 cwrk,cout,nfusemx,idoff,pspecl1,nchmx,nwmx,lpack,iwrk,ixs)
      print*,'CALL MKSEQ, NFUSEMX ',nfusemx
      call mkseq(ifdir,iset,nuse,nch,nd,
     1 cwrk,cout,nfusemx,idoff,pspecl1,nchmx,nwmx,lpack,iwrk,ixs,
     2 nsmax, x_scr, iuse_scr)
C     ME 	REMOVE SCRATCH FILE
c       close files; delete f_file on closing
      close(out_unit)
C ME 	REMOVE SCRATCH FILE
c      close(ioscr,status = 'delete')
C ME 	REMOVE SCRATCH FILE
500   continue

      if(.not.lseed) then
         print*,'FT another data file?'
         read(5,'(a1)') ans
         if(ans.eq.'y') go to 1000
      endif
      stop
      end

ext: f

max_stars_repo_name: simpeg-research/iris-mt-scratch

max_issues_repo_name: simpeg-research/iris-mt-scratch

max_forks_repo_name: simpeg-research/iris-mt-scratch


    subroutine red(value,hexrep,rfrac)
!-----------------------------------------------------------------------

! Provided with a real value ``value'' this subroutine returns the red
! portion of the color specification.

! Note common block "color" and call to "hexa"

!-----------------------------------------------------------------------
    real ::            value
    character       hexrep*2,hexa*2
    common /color/  cmin,cmax,cint(4),cscl
    hexrep='00'
    if(value < cint(1))then

    ! Scale it between (y0,0):

        integ=int((cint(1)-value)/(cint(1)-cmin)*cscl)
        if(integ > 255) integ = 255
        if(integ < 0)   integ = 0
    else if((value >= cint(1)) .AND. (value < cint(3)))then

    ! Scale it between (0,0):

        integ  = 0
    else if((value >= cint(3)) .AND. (value < cint(4)))then

    ! Scale it between (0,255):

        integ = int((value-cint(3))/(cint(4)-cint(3))*255.)
        if(integ > 255) integ = 255
        if(integ < 0)   integ = 0
    else if(value >= cint(4))then

    ! Scale it between (255,255):

        integ  = 255
    end if

! Establish coding and return:

    rfrac  = real(integ) / 255.
    hexrep = hexa(integ)
    return
    end subroutine red


text:
The Daily Democrat is the local Media newspaper in Woodland.

 Staff

Editor:
Jim Smith
news@dailydemocrat.com
(530) 4066230

Associate Editor:
Jake Dorsey
jdorsey@dailydemocrat.com
(530) 4066229

Reporter:
Crystal Lee
clee@dailydemocrat.com
(530) 4066232

Reporter:
Lizeth Cazares
lcazares@dailydemocrat.com
(530) 4066233

Letters to the Editor:
letters@dailydemocrat.com

Sports Editor:
Bruce Burton
sports@dailydemocrat.com
(530) 4066236

Photography:
photo@dailydemocrat.com

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
! This samples illustrates how to share scratch memory among plans to minimize memory usage.
!
! This can be done like this
! 1. Create (but don't make) plans
! 2. Call cufftSetAutoAllocation(plan, false) on all plans
! 3. Call cufftMakePlan3d(plan, ..., scratch_size) on all plans and retrieve the required scratch size per plan
! 4. Compute the maximum scratch size accros plans _AND_ accross MPI ranks (see note below on nvshmem_malloc)
! 5. Allocate memory using _cufftMpNvshmemMalloc
! 6. Call cufftSetWorkArea(plan, buffer) on all plans
! 7. Call cufftExec, cufftXtMemcpy, etc
! 8. Free memory using _cufftMpNvshmemFree
! 9. Destroy the plans

module cufft_required
    integer :: planr2c, planc2r
    integer :: local_rshape(3), local_rshape_permuted(3), local_permuted_cshape(3)

end module cufft_required

program cufftmp_r2c_workarea
    use iso_c_binding
    use cudafor
    use cufftXt
    use cufft
    use openacc
    use mpi
    use cufft_required
    implicit none

    integer :: size, rank, ndevices, ierr
    integer :: nx, ny, nz ! nx slowest
    integer :: i, j, k
    integer :: my_nx, my_ny, my_nz, ranks_cutoff, whichgpu(1)
    real, dimension(:, :, :), allocatable :: u, ref
    complex, dimension(:,:,:), allocatable :: u_permuted
    real :: max_norm, max_diff

    type(c_devptr)  :: workArea_ptr

    ! cufft stuff
    integer(c_size_t) :: worksize(2), localScratchSize, scratchSize
    type(cudaLibXtDesc), pointer :: u_desc
    type(cudaXtDesc), pointer    :: u_descptr
    complex, pointer, device     :: u_dptr(:,:,:)
    integer(kind=cuda_stream_kind) :: stream

    call mpi_init(ierr)
    call mpi_comm_size(MPI_COMM_WORLD,size,ierr)
    call mpi_comm_rank(MPI_COMM_WORLD,rank,ierr)

    call checkCuda(cudaGetDeviceCount(ndevices))
    call checkCuda(cudaSetDevice(mod(rank, ndevices)))
    whichgpu(1) = mod(rank, ndevices)

    print*,"Hello from rank ", rank, " gpu id", mod(rank, ndevices), "size", size

    nx = 256
    ny = nx
    nz = nx

    ! We start with X-Slabs
    ! Ranks 0 ... (nx % size - 1) have 1 more element in the X dimension
    ! and every rank own all elements in the Y and Z dimensions.
    ranks_cutoff = mod(nx, size)
    my_nx = nx / size
    if (rank < ranks_cutoff) my_nx = my_nx + 1
    my_ny =  ny;
    my_nz =  nz;
    local_rshape = [2*(nz/2+1), ny, my_nx]
    local_permuted_cshape = [nz/2+1, ny/size, nx]
    local_rshape_permuted = [2*(nz/2+1), ny/size, nx]
    if (mod(ny, size) > 0) then
        print*," ny has to divide evenly by mpi_procs"
        call mpi_finalize(ierr)
    end if
    if (rank == 0) then
        write(*,*) "local_rshape          :", local_rshape(1), local_rshape(2), local_rshape(3)
        write(*,*) "local_permuted_cshape :", local_permuted_cshape(1), local_permuted_cshape(2), local_permuted_cshape(3)
    end if

    ! Generate local, distributed data
    allocate(u(local_rshape(1), local_rshape(2), local_rshape(3)))
    allocate(u_permuted(local_permuted_cshape(1), local_permuted_cshape(2), local_permuted_cshape(3)))
    allocate(ref(local_rshape(1), local_rshape(2), local_rshape(3)))
    print*,'shape of u is ', shape(u)
    print*,'shape of u_permuted is ', shape(u_permuted)
    call generate_random(nz, local_rshape(1), local_rshape(2), local_rshape(3), u)
    ref = u
    u_permuted = (0.0,0.0)

    call checkNorm(nz, local_rshape(1), local_rshape(2), local_rshape(3), u, max_norm)
    print*, "initial data on ", rank, " max_norm is ", max_norm

    call checkCufft(cufftCreate(planr2c))
    call checkCufft(cufftCreate(planc2r))

#ifdef SHARED_WORKAREA
    call checkCufft(cufftSetAutoAllocation(planr2c, 0));
    call checkCufft(cufftSetAutoAllocation(planc2r, 0));
#endif
    call checkCufft(cufftMpAttachComm(planr2c, CUFFT_COMM_MPI, MPI_COMM_WORLD), 'cufftMpAttachComm error')
    call checkCufft(cufftMpAttachComm(planc2r, CUFFT_COMM_MPI, MPI_COMM_WORLD), 'cufftMpAttachComm error')

    call checkCufft(cufftMakePlan3d(planr2c, nz, ny, nx, CUFFT_R2C, worksize(0)), 'cufftMakePlan3d r2c error')
    call checkCufft(cufftMakePlan3d(planc2r, nz, ny, nx, CUFFT_C2R, worksize(1)), 'cufftMakePlan3d c2r error')

#ifdef SHARED_WORKAREA
    ! Compute how much scratch to allocate -- find the largest size for all the MPI ranks
    localScratchSize = max(worksize(0), worksize(1))
    write(*,'(A30,I2, A, F20.2,A)') 'cuFFT Plans allocated workarea for rank ', rank, ' : ', real(localScratchSize)/1024, " KB"
    call MPI_Allreduce(localScratchSize, scratchSize, 1, MPI_LONG_LONG_INT, MPI_MAX, MPI_COMM_WORLD, ierr)

    ! Allocate scratch size using NVSHMEM
    call checkCufft(cufftMpNvshmemMalloc(scratchSize, workArea_ptr))

    if (rank .eq. 0) then
        write(*,'(A30,F20.2,A)') 'cuFFT Plans allocated workarea on each GPU is : ', real(scratchSize)/1024, " KB"
    endif

    !  Pass the scratch to cuFFT
    call checkCufft(cufftSetWorkArea(planr2c, workArea_ptr))
    call checkCufft(cufftSetWorkArea(planc2r, workArea_ptr))
#endif
    call checkCufft(cufftXtMalloc(planr2c, u_desc, CUFFT_XT_FORMAT_INPLACE), 'cufftXtMalloc error')
    call checkCufft(cufftXtMemcpy(planr2c, u_desc, u, CUFFT_COPY_HOST_TO_DEVICE), 'cufftXtMemcpy error')
    !call cufft_memcpyH2D(u_desc, u, CUFFT_XT_FORMAT_INPLACE, .true.)
    ! now reset u to make sure the check later is valid
    u = 0.0

    !xxxxxxxxxxxxxxxxxxxxxxxxxx Forward
    !print*,2.5
    call checkCufft(cufftXtExecDescriptor(planr2c, u_desc, u_desc, CUFFT_FORWARD),'forward fft failed')
    ! in case we want to check the results after Forward
    !call checkCufft(cufftXtMemcpy(planr2c, u_permuted, u_desc, CUFFT_COPY_DEVICE_TO_HOST), 'permuted D2H error')
    !call checkNormComplex(local_permuted_cshape(1), local_permuted_cshape(2), local_permuted_cshape(3), u_permuted, max_norm)
    !write(*,'(A18, I1, A14, F25.8)') "after R2C ", rank, " max_norm is ", max_norm

    ! Data is now distributed as Y-Slab. We need to scale the output
    call c_f_pointer(u_desc%descriptor, u_descptr)

    ! print*,'doris rank ', rank, u_descptr%nGPUs, " " , u_descptr%GPUs(0)
    call c_f_pointer(u_descptr%data(1), u_dptr, [local_permuted_cshape(1), local_permuted_cshape(2), local_permuted_cshape(3)])
    !$cuf kernel do (3)
    do k =1, local_permuted_cshape(3)
        do j = 1, local_permuted_cshape(2)
            do i = 1, local_permuted_cshape(1)
                u_dptr(i,j,k) = u_dptr(i,j,k) / real(nx*ny*nz)
            end do
        end do
    end do
    call checkCuda(cudaDeviceSynchronize())

    ! in case we want to check again after scaling
    call checkCufft(cufftXtMemcpy(planr2c, u_permuted, u_desc, CUFFT_COPY_DEVICE_TO_HOST), 'permuted D2H error')
    call checkNormComplex(local_permuted_cshape(1), local_permuted_cshape(2), local_permuted_cshape(3), u_permuted, max_norm)
    write(*,'(A18, I1, A14, F25.8)') "after scaling ", rank, " max_norm is ", max_norm

    !xxxxxxxxxxxxxxxxxxxxxxxxxxxx inverse
    call checkCufft(cufftXtExecDescriptor(planc2r, u_desc, u_desc, CUFFT_INVERSE), 'inverse fft failed')
    call checkCufft(cufftXtMemcpy(planc2r, u, u_desc, CUFFT_COPY_DEVICE_TO_HOST), 'D2H failed')
    !call cufft_memcpyD2H(u, u_desc, CUFFT_XT_FORMAT_INPLACE, .true.)
    call checkCufft(cufftXtFree(u_desc))
#ifdef SHARED_WORKAREA
    call checkCufft(cufftMpNvshmemFree(workArea_ptr))
#endif
    call checkCufft(cufftDestroy(planr2c))
    call checkCufft(cufftDestroy(planc2r))

    call checkNormDiff(nz, local_rshape(1), local_rshape(2), local_rshape(3), u, ref, max_norm, max_diff)
    write(*,'(A18, I1, A14, F25.8, A14, F15.8)') "after C2R ", rank, " max_norm is ", max_norm, " max_diff is ", max_diff
    write(*,'(A25, I1, A14, F25.8)') "Relative Linf on rank ", rank, " is ", max_diff/max_norm
    deallocate(u)
    deallocate(ref)
    deallocate(u_permuted)

    call mpi_finalize(ierr)

    if(max_diff / max_norm > 1e-5) then
        print*, ">>>> FAILED on rank ", rank
        stop (1)
    else
        print*, ">>>> PASSED on rank ", rank
    end if

contains
    subroutine checkCuda(istat, message)
        implicit none
        integer, intent(in)                   :: istat
        character(len=*),intent(in), optional :: message
        if (istat /= cudaSuccess) then
            write(*,"('Error code: ',I0, ': ')") istat
            write(*,*) cudaGetErrorString(istat)
            if(present(message)) write(*,*) message
            call mpi_finalize(ierr)
        endif
    end subroutine checkCuda

    subroutine checkCufft(istat, message)
        implicit none
        integer, intent(in)                   :: istat
        character(len=*),intent(in), optional :: message
        if (istat /= CUFFT_SUCCESS) then
            write(*,"('Error code: ',I0, ': ')") istat
            write(*,*) cudaGetErrorString(istat)
            if(present(message)) write(*,*) message
            call mpi_finalize(ierr)
        endif
    end subroutine checkCufft

    subroutine generate_random(nz1, nz, ny, nx, data)
        implicit none
        integer, intent(in) :: nx, ny, nz, nz1
        real, dimension(nz, ny, nx), intent(out) :: data
        real :: rand(1)
        integer :: i,j,k
        !call random_seed(put=(/seed, seed+1/))
        do k =1, nx
            do j = 1, ny
                do i = 1, nz1
                    call random_number(rand)
                    data(i,j,k) = rand(1)
                end do
            end do
        end do

    end subroutine generate_random

    subroutine checkNorm(nz1, nz, ny, nx, data, max_norm)
        implicit none
        integer, intent(in)  :: nx, ny, nz, nz1
        real, dimension(nz, ny, nx), intent(in) :: data
        real :: max_norm
        integer :: i, j, k
        max_norm = 0
        do k =1, nx
            do j = 1, ny
                do i = 1, nz1
                    max_norm = max(max_norm, abs(data(i,j,k)))
                    !write(*,'(A10 I3 I3 I3 A2 F18.8 A10 I2)') "ref[", &
                    !i, j, k, "]", ref(i,j,k), " at rank ", rank
                end do
            end do
        end do
    end subroutine checkNorm

    subroutine checkNormComplex(nz, ny, nx, data, max_norm)
        implicit none
        integer, intent(in)  :: nx, ny, nz
        complex, dimension(nz, ny, nx), intent(in) :: data
        real :: max_norm, max_diff
        integer :: i,j,k
        max_norm = 0
        do k =1, nx
            do j = 1, ny
                do i = 1, nz
                    max_norm = max(max_norm, abs(data(i,j,k)%re))
                    max_norm = max(max_norm, abs(data(i,j,k)%im))
                end do
            end do
        end do
    end subroutine checkNormComplex

    subroutine checkNormDiff(nz1, nz, ny, nx, data, ref, max_norm, max_diff)
        implicit none
        integer, intent(in)  :: nx, ny, nz, nz1
        real, dimension(nz, ny, nx), intent(in) :: data, ref
        real :: max_norm, max_diff
        max_norm = 0
        max_diff = 0
        do k =1, nx
            do j = 1, ny
                do i = 1, nz1
                    max_norm = max(max_norm, abs(data(i,j,k)))
                    max_diff = max(max_diff, abs(ref(i,j,k)-data(i,j,k)))
                    !if (abs(ref(i,j,k)-data(i,j,k)) > 0.0001) write(*,'(A10 I3 I3 I3 A2 F18.8 A7 F18.8 A10 I2)') "diff ref[", &
                    !i, j, k, "]", ref(i,j,k), "data ", data(i,j,k), " at rank ", rank
                end do
            end do
        end do
    end subroutine checkNormDiff

#ifdef ACC
    subroutine scalingData(nz, ny, nx, data, factor)
        implicit none
        integer, intent(in)  :: nx, ny, nz
        complex, dimension(nz, ny, nz) :: data
        !$acc declare deviceptr(data)
        real, intent(in)     :: factor

    !$acc parallel loop collapse(3)
        do k =1, nx
            do j = 1, ny
                do i = 1, nz
                    data(i, j, k) = data(i, j, k) / factor
                end do
            end do
        end do

    end subroutine scalingData
#endif

subroutine cufft_memcpyH2D(ulibxt, u_h, data_format, ismemcpy)
    implicit none
    type(cudaLibXtDesc), pointer, intent(out) :: ulibxt
    real, dimension(*), intent(in)          :: u_h
    integer, intent(in)                         :: data_format
    logical, intent(in)                         :: ismemcpy
    type(cudaXtDesc), pointer  :: uxt
    real, dimension(:,:,:), device, pointer :: u_d

    if(data_format == CUFFT_XT_FORMAT_INPLACE_SHUFFLED) then
      if (ismemcpy == .false.) then
        call checkCufft(cufftXtMemcpy(planc2r, ulibxt, u_h, CUFFT_COPY_HOST_TO_DEVICE), "cufft_memcpyHToD pinv Error")
      else
        call c_f_pointer(ulibxt%descriptor, uxt)
        call c_f_pointer(uxt%data(1), u_d, local_rshape_permuted)
        call checkCuda(cudaMemcpy(u_d, u_h, product(int(local_rshape_permuted,kind=8))), "cudamemcpy H2D Error")
        nullify(u_d, uxt)
      endif
    endif

    if (data_format == CUFFT_XT_FORMAT_INPLACE) then
      if (ismemcpy == .false.) then
        call checkCufft(cufftXtMemcpy(planr2c, ulibxt, u_h, CUFFT_COPY_HOST_TO_DEVICE), "cufft_memcpyHToD pfor Error")
      else
        call c_f_pointer(ulibxt%descriptor, uxt)
        call c_f_pointer(uxt%data(1), u_d, local_rshape)
        call checkCuda(cudaMemcpy(u_d, u_h, product(int(local_rshape,kind=8))), "cudamemcpy H2D Error")
        nullify(u_d, uxt)
      endif
    endif
end subroutine cufft_memcpyH2D

subroutine cufft_memcpyD2H(u_h, ulibxt, data_format,ismemcpy)
    implicit none
    type(cudaLibXtDesc), pointer, intent(in) :: ulibxt
    real, dimension(*), intent(out)      :: u_h
    integer, intent(in)                      :: data_format
    logical, intent(in)                      :: ismemcpy
    type(cudaXtDesc), pointer  :: uxt
    real, dimension(:,:,:), device, pointer :: u_d

    if(data_format == CUFFT_XT_FORMAT_INPLACE_SHUFFLED) then
      if (ismemcpy == .false.) then
        call checkCufft(cufftXtMemcpy(planr2c, u_h, ulibxt, CUFFT_COPY_DEVICE_TO_HOST), "cufft_memcpyDToH pfor Error")
      else
        call c_f_pointer(ulibxt%descriptor, uxt)
        call c_f_pointer(uxt%data(1), u_d, local_rshape_permuted)
        call checkCuda(cudaMemcpy(u, u_d, product(int(local_rshape_permuted,kind=8))), "cudamemcpy D2H Error")
        nullify(u_d, uxt)
      endif
    endif

    if (data_format == CUFFT_XT_FORMAT_INPLACE) then
      if (ismemcpy == .false.) then
        call checkCufft(cufftXtMemcpy(planc2r, u_h, ulibxt, CUFFT_COPY_DEVICE_TO_HOST), "cufft_memcpyDToH pinv Error")
      else
        call c_f_pointer(ulibxt%descriptor, uxt)
        call c_f_pointer(uxt%data(1), u_d, local_rshape)
        call checkCufft(cudamemcpy(u, u_d, product(int(local_rshape,kind=8))), "cufft_memcpyD2H error")
        nullify(u_d, uxt)
      endif
    endif
end subroutine cufft_memcpyD2H

end program cufftmp_r2c_workarea

max_stars_repo_name: katherineding/CUDALibrarySamples

max_issues_repo_name: katherineding/CUDALibrarySamples

max_forks_repo_name: katherineding/CUDALibrarySamples


!
! drawButton.f90
! Fortran から　libui　を使う。
!
! Author Yoji.Hosokawa 2016/4
!
! gfortran -O2 -g -Wall -o ./drawLabel_center drawLabel_center.f90 -L. -lui -llibuif_parts `pkg-config gtk+-3.0 --libs` -ldl -lpthread
!

module handlers
  use iso_c_binding
  use libuif_structs
  use libuif
  implicit none

    type(c_ptr) :: drawLabel,drawButton
    type(uiAreaHandler) :: drawLabelHandler,drawButtonHandler

    character(200) :: drawLabelText,drawButtonText
    integer :: onButton = 0, onClick = 0

!// and some colors
!// names and values from https://msdn.microsoft.com/en-us/library/windows/desktop/dd370907%28v=vs.85%29.aspx
    integer,parameter :: colorWhite = x'FFFFFF'
    integer,parameter :: colorBlack = x'000000'
    integer,parameter :: colorRed = x'FF0000'
    integer,parameter :: colorDodgerBlue = x'1E90FF'
    integer,parameter :: d2dBlack = x'000000'
    integer,parameter :: d2dLightSlateGray = x'778899'
    integer,parameter :: d2dCornflowerBlue = x'6495ED'
    integer,parameter :: d2dWhite = x'FFFFFF'
    integer,parameter :: d2dYellowGreen = x'9ACD32'
    integer,parameter :: d2dYellow = x'FFFF00'
    integer,parameter :: d2dForestGreen = x'228B22'
    integer,parameter :: d2dOliveDrab = x'6B8E23'
    integer,parameter :: d2dLightSkyBlue = x'87CEFA'

contains

! "close mainwindow" is a this app
  subroutine onClosing()
    use iso_c_binding
    use libuif
      call uiQuit()
  end subroutine onClosing

  subroutine drawLabelup()
    use iso_c_binding
    use libuif
      drawLabelText = "押されました！！"//nilc
      call uiAreaQueueRedrawAll(drawLabel)
  end subroutine drawLabelup

  subroutine d2dColorToRGB(color, r, g, b)
    use iso_c_binding
    use libuif
      integer :: color
      double precision :: r,g,b
        r = dble(ibits(color,16,8)) / 255.0d0
        g = dble(ibits(color,8,8)) / 255.0d0
        b = dble(ibits(color,0,8)) / 255.0d0
  end subroutine d2dColorToRGB

  subroutine d2dSolidBrush(brush, color, alpha)
    use iso_c_binding
    use libuif_structs
    use libuif
      type(uiDrawBrush) :: brush
      integer :: color
      double precision :: alpha
        brush%Type = uiDrawBrushTypeSolid
        call d2dColorToRGB(color, brush%R, brush%G, brush%B)
        brush%A = alpha;
  end subroutine d2dSolidBrush

  subroutine handlerDraw(a, area, dp)
    use iso_c_binding
    use libuif_structs
    use libuif
    use libuif_parts
      implicit none
      type(uiAreaHandler) :: a
      type(c_ptr) :: area
      type(c_ptr) :: path
      type(uiDrawBrush) :: brush
      type(uiAreaDrawParams) :: dp
      type(c_ptr) :: layout
      type(c_ptr) :: textFont

      double precision :: cx,cy,fontSize,strlenw
      type(uiDrawTextFontDescriptor),pointer :: font
      character(100,c_char),pointer :: family(:)
      allocate(family(1))
      allocate(font)

!  fill the area with white
       call d2dSolidBrush(brush, colorWhite, 1.0d0)
       path = uiDrawNewPath(uiDrawFillModeWinding)
       call uiDrawPathAddRectangle(path, 0.0_c_double, 0.0_c_double, dp%AreaWidth, dp%AreaHeight)
       call uiDrawPathEnd(path)
       call uiDrawFill(dp%Context, path, brush)
       call uiDrawFreePath(path)

! dummy draw with color in area and set foreground color
       call d2dSolidBrush(brush, colorBlack, 1.0d0)
       path = uiDrawNewPath(uiDrawFillModeAlternate)
       call uiDrawPathEnd(path)
       call uiDrawFill(dp%Context, path, brush)
       call uiDrawFreePath(path)

! select font size & position
      fontSize = 20.0d0
      strlenw = utf8lenw(drawLabelText)

      cx = dp%AreaWidth / 2.0d0 - fontSize * strlenw * 0.64d0 / 2.0d0
      if (cx < 0) then
        cx = 0.0d0
      end if
      cy = dp%AreaHeight / 2.0d0 - fontSize / 2.0d0 - fontSize * 0.15d0

! set font
!      family(1) = "Monospace"//nilc
      family(1) = "Arial"//nilc
      font%Family = c_loc(family(1))
      font%Size = fontSize
      font%Weight = uiDrawTextWeightNormal
      font%Italic = uiDrawTextItalicNormal
      font%Stretch = uiDrawTextStretchNormal
      textFont = uiDrawLoadClosestFont(c_loc(font))

! draw text

      layout = uiDrawNewTextLayout(drawLabelText, textFont, -1.0d0)
      call uiDrawFreeTextFont(textFont)

      call uiDrawText(dp%Context, cx, cy , layout)
      call uiDrawFreeTextLayout(layout)

  end subroutine handlerDraw

  subroutine handlerMouseEvent()
!  do nothing
  end subroutine handlerMouseEvent

  subroutine handlerMouseCrossed()
    use iso_c_binding
    use libuif
!  do nothing
  end subroutine handlerMouseCrossed

  subroutine handlerDragBroken()
    use iso_c_binding
    use libuif
!  do nothing
  end subroutine handlerDragBroken

  function handlerKeyEvent() result(ret)
    use iso_c_binding
    use libuif
      integer :: ret
!  reject all keys
    ret = 0
  end function handlerKeyEvent

  subroutine handlerDrawButton(a, area, dp)
    use iso_c_binding
    use libuif_structs
    use libuif
    use libuif_parts
      implicit none
      type(uiAreaHandler) :: a
      type(c_ptr) :: area
      type(c_ptr) :: path
      type(uiDrawBrush) :: brush
      type(uiAreaDrawParams) :: dp
      type(c_ptr) :: layout
      type(c_ptr) :: textFont

      double precision :: cx,cy,fontSize,strlenw
      type(uiDrawTextFontDescriptor),pointer :: font
      character(100,c_char),pointer :: family(:)
      allocate(family(1))
      allocate(font)

!  fill the area with d2dLightSlateGray
       if ( onButton == 1) then
         call d2dSolidBrush(brush, d2dLightSlateGray, 0.3d0)
         if (onClick == 1) then
           call d2dSolidBrush(brush, d2dLightSlateGray, 0.7d0)
         else
           call d2dSolidBrush(brush, d2dLightSlateGray, 0.3d0)
         endif
       else
         call d2dSolidBrush(brush, d2dLightSlateGray, 0.1d0)
       endif

       path = uiDrawNewPath(uiDrawFillModeWinding)
       call uiDrawPathAddRectangle(path, 0.0_c_double, 0.0_c_double, dp%AreaWidth, dp%AreaHeight)
       call uiDrawPathEnd(path)
       call uiDrawFill(dp%Context, path, brush)
       call uiDrawFreePath(path)

! dummy draw with color in area and set foreground color
       call d2dSolidBrush(brush, colorBlack, 1.0d0)
       path = uiDrawNewPath(uiDrawFillModeAlternate)
       call uiDrawPathEnd(path)
       call uiDrawFill(dp%Context, path, brush)
       call uiDrawFreePath(path)

! select font size & position
      fontSize = 20.0d0
      strlenw = utf8lenw(drawButtonText)

      cx = dp%AreaWidth / 2.0d0 - fontSize * strlenw * 0.64d0 / 2.0d0
      if (cx < 0) then
        cx = 0.0d0
      end if
      cy = dp%AreaHeight / 2.0d0 - fontSize / 2.0d0 - fontSize * 0.15d0

! set font
!      family(1) = "Monospace"//nilc
      family(1) = "Arial"//nilc
      font%Family = c_loc(family(1))
      font%Size = fontSize
      font%Weight = uiDrawTextWeightNormal
      font%Italic = uiDrawTextItalicNormal
      font%Stretch = uiDrawTextStretchNormal
      textFont = uiDrawLoadClosestFont(c_loc(font))

! draw text

      layout = uiDrawNewTextLayout(drawButtonText, textFont, -1.0d0)
      call uiDrawFreeTextFont(textFont)

      call uiDrawText(dp%Context, cx, cy , layout)
      call uiDrawFreeTextLayout(layout)

  end subroutine handlerDrawButton

  subroutine drawButtonMouseEvent(a, area, e)
    use iso_c_binding
    use libuif_structs
    use libuif
    implicit none
      type(uiAreaHandler) :: a
      type(c_ptr) :: area
      type(uiAreaMouseEvent) :: e

      if (e%X >= 5.0d0 .and. e%X <= e%AreaWidth - 5.0d0 .and. e%Y >= 5.0d0 .and. e%Y <= e%AreaHeight - 5.0d0) then
        onButton = 1
        call uiAreaQueueRedrawAll(drawButton)
      else
        onButton = 0
        call uiAreaQueueRedrawAll(drawButton)
      endif

! when mouse clicked
! left: 1, center: 2, right: 3
      if ( e%Down == 1) then
        onClick = 1
        call uiAreaQueueRedrawAll(drawButton)
        call drawLabelup()
      else
        onClick = 0
        call uiAreaQueueRedrawAll(drawButton)
      end if

  end subroutine drawButtonMouseEvent

  subroutine drawButtonMouseCrossed()
    use iso_c_binding
    use libuif
!  do nothing
  end subroutine drawButtonMouseCrossed

  subroutine drawButtonDragBroken()
    use iso_c_binding
    use libuif
!  do nothing
  end subroutine drawButtonDragBroken

  function drawButtonKeyEvent() result(ret)
    use iso_c_binding
    use libuif
      integer :: ret
!  reject all keys
    ret = 0
  end function drawButtonKeyEvent

end module handlers

program main
    use iso_c_binding
    use libuif_structs
    use libuif
    use handlers

    implicit none
    type(c_ptr) :: mainwin
    type(c_ptr) :: vbox,ebtn
    type(c_ptr) :: msggrid,hbox

! initialize
    call uiInit(0)

! start procedure

    mainwin = uiNewWindow("drawLabel centering Test"//c_null_char, 640, 480, 0)
    call uiWindowSetMargined(mainwin, 1)
    call uiWindowOnClosing(mainwin, c_funloc(onClosing), c_null_ptr)

    vbox = uiNewVerticalBox()
    call uiWindowSetChild(mainwin, vbox)

    hbox = uiNewHorizontalBox()
    call uiBoxAppend(vbox, hbox, 1)

    msggrid = uiNewGrid()
    call uiGridSetPadded(msggrid, 1)
    call uiBoxAppend(hbox, msggrid, 1)

! drawLabel
    drawLabelHandler%Draw = c_funloc(handlerDraw)
    drawLabelHandler%MouseEvent = c_funloc(handlerMouseEvent)
    drawLabelHandler%MouseCrossed = c_funloc(handlerMouseCrossed)
    drawLabelHandler%DragBroken = c_funloc(handlerDragBroken)
    drawLabelHandler%KeyEvent = c_funloc(handlerKeyEvent)
    drawLabel = uiNewArea(drawLabelHandler)
    call uiGridAppend(msggrid, drawLabel,0,0,2,1,1, uiAlignFill, 1, uiAlignFill)
    drawLabelText = "真ん中に表示するかな？？？"//nilc

! drawButton
    drawButtonHandler%Draw = c_funloc(handlerDrawButton)
    drawButtonHandler%MouseEvent = c_funloc(drawButtonMouseEvent)
    drawButtonHandler%MouseCrossed = c_funloc(drawButtonMouseCrossed)
    drawButtonHandler%DragBroken = c_funloc(drawButtonDragBroken)
    drawButtonHandler%KeyEvent = c_funloc(drawButtonKeyEvent)
    drawButton = uiNewArea(drawButtonHandler)
    call uiGridAppend(msggrid, drawButton,0,2,2,1,1, uiAlignFill, 1, uiAlignFill)
    drawButtonText = "ドローボタンです。イベントは使えるか？"//nilc

!    call uiBoxSetPadded(vbox,1)
    ebtn = uiNewButton("終了"//c_null_char)
    call uiButtonOnClicked(ebtn, c_funloc(onClosing), c_null_ptr)
    call uiGridAppend(msggrid, ebtn,1,4,1,1,1, uiAlignFill, 0, uiAlignFill)

    call uiControlShow(mainwin)
    call uiMain()

end program main


text:
        subroutine idz_estrank0(eps,m,n,a,w,n2,krank,ra,rat,scal)
c
c       routine idz_estrank serves as a memory wrapper
c       for the present routine. (Please see routine idz_estrank
c       for further documentation.)
c
        implicit none
        integer m,n,n2,krank,ifrescal,k,nulls
        real*8 eps,scal(n2+1)
        complex*16 a(m,n),ra(n2,n),residual,w(17*m+70),rat(n,n2+1)
c
c
c       Apply the random matrix to every column of a, obtaining ra.
c
        do k = 1,n
          call idz_frm(m,n2,w,a(1,k),ra(1,k))
        enddo ! k
c
c
c       Transpose ra to obtain rat.
c
        call idz_transposer(n2,n,ra,rat)
c
c
        krank = 0
        nulls = 0
c
c
c       Loop until nulls = 7, krank+nulls = n2, or krank+nulls = n.
c
 1000   continue
c
c
          if(krank .gt. 0) then
c
c           Apply the previous Householder transformations
c           to rat(:,krank+1).
c
            ifrescal = 0
c
            do k = 1,krank
              call idz_houseapp(n-k+1,rat(1,k),rat(k,krank+1),
     1                          ifrescal,scal(k),rat(k,krank+1))
            enddo ! k
c
          endif ! krank .gt. 0
c
c
c         Compute the Householder vector associated
c         with rat(krank+1:*,krank+1).
c
          call idz_house(n-krank,rat(krank+1,krank+1),
     1                   residual,rat(1,krank+1),scal(krank+1))
c
c
          krank = krank+1
          if(abs(residual) .le. eps) nulls = nulls+1
c
c
        if(nulls .lt. 7 .and. krank+nulls .lt. n2
     1   .and. krank+nulls .lt. n)
     2   goto 1000
c
c
        if(nulls .lt. 7) krank = 0
c
c
        return
        end
c
c
c
c


text:
      program rsac
      implicit none

!     Define the Maximum size of the data Array
      integer MAX
      parameter (MAX=1000)

!     Define the Data Array of size MAX
      real yarray
      dimension yarray(MAX)

!     Declare Variables used in the rsac1() subroutine
      real beg, del
      integer nlen
      character*10 KNAME
      integer nerr

!     Define the file to be read      g
      kname = 'FILE1'

!     Call rsac1 to read filename kname
!        - Data is loaded into yarray
!        - Length of data is stored in nlen
!        - Begining time and time sampling are in beg and del
!        - MAX is the maximum number of points to be read in
!        - nerr is the Error return flag
      call rsac1(kname, yarray, nlen, beg, del, MAX, nerr)

!     Check the error status, nerr
!        - 0 on Success
!        - Non-Zero on Failure
      if(nerr .NE. 0) then
      	  write(*,*)'Error reading in file: ',kname
	  call exit(-1)
      endif

!     Do some processing ....

      call exit(0)
      end


text:
      SUBROUTINE RADF3 (IDO,L1,CC,CH,WA1,WA2)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION       CH(IDO,3,L1)           ,CC(IDO,L1,3)           ,
     1                WA1(*)     ,WA2(*)
C     *** TAUI IS -SQRT(3)/2 ***
      DATA TAUR,TAUI /-0.5D0,0.86602540378443864676D0/
      DO 101 K=1,L1
         CR2 = CC(1,K,2)+CC(1,K,3)
         CH(1,1,K) = CC(1,K,1)+CR2
         CH(1,3,K) = TAUI*(CC(1,K,3)-CC(1,K,2))
         CH(IDO,2,K) = CC(1,K,1)+TAUR*CR2
  101 CONTINUE
      IF (IDO .EQ. 1) RETURN
      IDP2 = IDO+2
      DO 103 K=1,L1
         DO 102 I=3,IDO,2
            IC = IDP2-I
            DR2 = WA1(I-2)*CC(I-1,K,2)+WA1(I-1)*CC(I,K,2)
            DI2 = WA1(I-2)*CC(I,K,2)-WA1(I-1)*CC(I-1,K,2)
            DR3 = WA2(I-2)*CC(I-1,K,3)+WA2(I-1)*CC(I,K,3)
            DI3 = WA2(I-2)*CC(I,K,3)-WA2(I-1)*CC(I-1,K,3)
            CR2 = DR2+DR3
            CI2 = DI2+DI3
            CH(I-1,1,K) = CC(I-1,K,1)+CR2
            CH(I,1,K) = CC(I,K,1)+CI2
            TR2 = CC(I-1,K,1)+TAUR*CR2
            TI2 = CC(I,K,1)+TAUR*CI2
            TR3 = TAUI*(DI2-DI3)
            TI3 = TAUI*(DR3-DR2)
            CH(I-1,3,K) = TR2+TR3
            CH(IC-1,2,K) = TR2-TR3
            CH(I,3,K) = TI2+TI3
            CH(IC,2,K) = TI3-TI2
  102    CONTINUE
  103 CONTINUE
      RETURN
      END

max_stars_repo_name: utastudents/selalib

max_issues_repo_name: utastudents/selalib

max_forks_repo_name: utastudents/selalib