text:
subroutine z_difuflux(stage   ,lundia    ,kmax      ,nmmax     ,nmmaxj    , &
                  & lstsci    ,r0        ,r1        ,qxk       ,qyk       , &
                  & u         ,v         , &
                  & dicuv     ,guv       ,gvu       ,areau     ,areav     , &
                  & kfuz1     ,kfvz1     ,kfsz1     ,kcs       ,kfs       , &
                  & kfu       ,kfuz0     ,kfv       ,kfvz0     , &
                  & kfsmx0    ,kfsmax    ,kfsz0     , &
                  & kfumin    ,kfumx0    ,kfvmin    ,kfvmx0    ,sigdif    , &
                  & timest    ,icx       ,icy       ,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_difuflux.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_difuflux.f90 $
!!--description-----------------------------------------------------------------
!
!    Function: Compute flux corresponding to Z_DIFU
!
! Method used:
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use mathconsts
    !
    use globaldata
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    integer                             , pointer :: lsed
    real(fp) , dimension(:,:,:)         , pointer :: fluxu
    real(fp) , dimension(:,:,:)         , pointer :: fluxuc
    real(fp) , dimension(:,:,:)         , pointer :: fluxv
    real(fp) , dimension(:,:,:)         , pointer :: fluxvc
    character(13)                       , pointer :: trasol
    type (flwoutputtype)                , pointer :: flwoutput
    type (gd_flwpar)                    , pointer :: gdflwpar
!
! Global variables
!
    integer                                                       , intent(in) :: icx    !  Description and declaration in inout.igs
    integer                                                       , intent(in) :: icy    !  Description and declaration in inout.igs
    integer                                                       , intent(in) :: kmax   !  Description and declaration in inout.igs
    integer                                                       , intent(in) :: nmmax  !  Description and declaration in inout.igs
    integer                                                       , intent(in) :: nmmaxj !  Description and declaration in inout.igs
    integer                                                       , intent(in) :: lstsci !  Description and declaration in inout.igs
    integer                                                       , intent(in) :: lundia !  Description and declaration in inout.igs
    integer      , dimension(gdp%d%nmlb:gdp%d%nmub)               , intent(in) :: kcs    !  Description and declaration in esm_alloc_int.f90

    integer      , dimension(gdp%d%nmlb:gdp%d%nmub)               , intent(in) :: kfu    !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(gdp%d%nmlb:gdp%d%nmub)               , intent(in) :: kfumin !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(gdp%d%nmlb:gdp%d%nmub)               , intent(in) :: kfumx0 !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)         , intent(in) :: kfuz0  !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(gdp%d%nmlb:gdp%d%nmub)               , intent(in) :: kfv    !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(gdp%d%nmlb:gdp%d%nmub)               , intent(in) :: kfvmin !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(gdp%d%nmlb:gdp%d%nmub)               , intent(in) :: kfvmx0 !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)         , intent(in) :: kfvz0  !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(gdp%d%nmlb:gdp%d%nmub)               , intent(in) :: kfs    !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(gdp%d%nmlb:gdp%d%nmub)               , intent(in) :: kfsmx0
    integer      , dimension(gdp%d%nmlb:gdp%d%nmub)               , intent(in) :: kfsmax
    integer      , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)         , intent(in) :: kfsz0  !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)         , intent(in) :: kfsz1  !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)         , intent(in) :: kfuz1  !  Description and declaration in esm_alloc_int.f90
    integer      , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)         , intent(in) :: kfvz1  !  Description and declaration in esm_alloc_int.f90
    real(fp)                                                      , intent(in) :: timest !!  Half Integration time step [sec.]
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)         , intent(in) :: areau
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)         , intent(in) :: areav
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)         , intent(in) :: dicuv  !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)               , intent(in) :: guv    !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub)               , intent(in) :: gvu    !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)         , intent(in) :: qxk    !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)         , intent(in) :: qyk    !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax, lstsci) , intent(in) :: r0     !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax, lstsci) , intent(in) :: r1     !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(lstsci)                              , intent(in) :: sigdif !  Description and declaration in esm_alloc_real.f90
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)         , intent(in) :: u
    real(fp)     , dimension(gdp%d%nmlb:gdp%d%nmub, kmax)         , intent(in) :: v
    character(8)                                                  , intent(in) :: stage  !!  First or second half time step
!
! Local variables
!
    integer  :: istat
    integer  :: k
    integer  :: l
    integer  :: kenu, kenv
    integer  :: kmin
    integer  :: kkmin, kkminu, kkminv
    integer  :: mink2
    integer  :: mink, mink_old, mink_new
    integer  :: ndm
    integer  :: nm
    integer  :: nmd
    integer  :: nmu
    integer  :: nmuu
    integer  :: num
    integer  :: nuum
    real(fp) :: cl, cr, difl, difr
    real(fp) :: cfl
    real(fp) :: flux
    real(fp) :: r00
    real(fp) :: rr1, rr2
    real(fp) :: rmax
    real(fp) :: rmin
    real(fp) :: switch

!
!! executable statements -------------------------------------------------------
!
    lsed           => gdp%d%lsed
    fluxu          => gdp%gdflwpar%fluxu
    fluxuc         => gdp%gdflwpar%fluxuc
    fluxv          => gdp%gdflwpar%fluxv
    fluxvc         => gdp%gdflwpar%fluxvc
    trasol         => gdp%gdtricom%trasol
    flwoutput      => gdp%gdflwpar%flwoutput
    gdflwpar       => gdp%gdflwpar

    !
    if (.not. flwoutput%difuflux .and. lsed == 0) return
    !
    istat = 0
    if (.not. associated(gdflwpar%fluxu)) then
       if (istat==0) allocate (gdflwpar%fluxu(gdp%d%nmlb:gdp%d%nmub, kmax, lstsci), stat = istat)
       if (istat==0) allocate (gdflwpar%fluxv(gdp%d%nmlb:gdp%d%nmub, kmax, lstsci), stat = istat)
       if (flwoutput%cumdifuflux) then
          if (istat==0) allocate (gdflwpar%fluxuc(gdp%d%nmlb:gdp%d%nmub, kmax, lstsci), stat = istat)
          if (istat==0) allocate (gdflwpar%fluxvc(gdp%d%nmlb:gdp%d%nmub, kmax, lstsci), stat = istat)
       endif
       !
       if (istat /= 0) then
          call prterr(lundia, 'U021', 'Z_DIFUFLUX: memory alloc error')
          call d3stop(1, gdp)
       endif
       !
       ! Define GDP pointers after the GDP allocations
       !
       fluxu          => gdp%gdflwpar%fluxu
       fluxv          => gdp%gdflwpar%fluxv
       flwoutput      => gdp%gdflwpar%flwoutput
       gdflwpar       => gdp%gdflwpar
       !
       if (flwoutput%cumdifuflux) then
          fluxuc => gdp%gdflwpar%fluxuc
          fluxvc => gdp%gdflwpar%fluxvc
          fluxuc = 0.0_fp
          fluxvc = 0.0_fp
       endif
    endif
    !
    ! Initialization
    !
    fluxu = 0.0_fp
    fluxv = 0.0_fp
    !

    if (trasol == 'van leer-2   ') then

       do nm = 1, nmmax
          nmd  = nm  - icx
          nmu  = nm  + icx
          nmuu = nmu + icx
          ndm  = nm  - icy
          num  = nm  + icy
          nuum = num + icy
          !
          ! ADVECTIVE TRANSPORT IN X-DIRECTION
          !
          if (kfu(nm) /= 0) then
             kmin = max(kfumin(nm), 1)
             do k = kmin, kmax
                if (kfuz0(nm,k) == 1) then
                   cfl = u(nm, k) * timest / gvu(nm)
                   if (qxk(nm, k) > 0.0_fp) then
                      do l = 1, lstsci
                         rr1 = abs(r0(nmd, k, l) - 2.0_fp*r0(nm, k, l) + r0(nmu, k, l))
                         rr2 = abs(r0(nmd, k, l) - r0(nmu, k, l))
                         if (kfsz0(nmd, k)*kfsz0(nmu, k) == 0 .or. rr1 >= rr2 .or. rr2 < eps_fp .or. kcs(nm)==3 .or. kcs(nmu)==3) then
                            r00 = r0(nm, k, l)
                         else
                            r00 = r0(nm , k, l)                                     &
                                & + (1.0_fp - cfl)*(r0(nm , k, l) - r0(nmd, k, l))  &
                                &                 *(r0(nmu, k, l) - r0(nm , k, l))  &
                                &                 /(r0(nmu, k, l) - r0(nmd, k, l))
                         endif
                         flux = qxk(nm, k) * r00
                         if (kcs(nm) == 1) then
                            fluxu(nm, k, l) = fluxu(nm, k, l) + flux
                         endif
                      enddo
                   elseif (qxk(nm, k) < 0.0_fp) then
                      do l = 1, lstsci
                         rr1 = abs(r0(nmuu, k, l) - 2.0_fp*r0(nmu, k, l) + r0(nm, k, l))
                         rr2 = abs(r0(nmuu, k, l) - r0(nm, k, l))
                         if (kfsz0(nm, k)*kfsz0(nmuu, k) == 0 .or. rr1 >= rr2 .or. rr2 < eps_fp .or. kcs(nm)==3 .or. kcs(nmu)==3) then
                            r00 = r0(nmu, k, l)
                         else
                            r00 = r0(nmu, k, l)                                      &
                                & + (1.0_fp + cfl)*(r0(nm , k, l) - r0(nmu , k, l))  &
                                &                 *(r0(nmu, k, l) - r0(nmuu, k, l))  &
                                &                 /(r0(nm , k, l) - r0(nmuu, k, l))
                         endif
                         flux = qxk(nm, k) * r00
                         if (kcs(nm) == 1) then
                            fluxu(nm, k, l) = fluxu(nm, k, l) + flux
                         endif
                      enddo
                   endif
                endif
             enddo
          endif
          !
          ! DIFFUSIVE TRANSPORT IN X-DIRECTION
          !
          if (kfs(nm)*kfs(nmu) == 1) then
             do k = kfumin(nm), kmax
                if (kfuz0(nm, k)==1 .and. kfsz0(nm, k)*kfsz0(nmu, k)==1) then
                   difl = dicuv(nm, k)
                   difr = dicuv(nmu, k)
                   kenu = max(0, abs(2 - kcs(nm)))
                   do l = 1, lstsci
                      cl              = r0(nm, k, l)
                      cr              = r0(nmu, k, l)
                      flux            = 0.5_fp*(cr - cl)*(difl + difr)/(sigdif(l)*gvu(nm))
                      fluxu(nm, k, l) = fluxu(nm, k, l) - areau(nm, k)*flux*kenu
                   enddo
                endif
             enddo
          endif
          !
          ! ADVECTIVE TRANSPORT IN Y-DIRECTION
          !
          if (kfv(nm) /= 0) then
             kmin = max(kfvmin(nm), 1)
             do k = kmin, kmax
                if (kfvz0(nm,k) == 1) then
                   cfl = v(nm, k) * timest / guv(nm)
                   if (qyk(nm, k) > 0.0_fp) then
                      do l = 1, lstsci
                         rr1 = abs(r0(ndm, k, l) - 2.0_fp*r0(nm, k, l) + r0(num, k, l))
                         rr2 = abs(r0(ndm, k, l) - r0(num, k, l))
                         if (kfsz0(ndm, k)*kfsz0(num, k) == 0 .or. rr1 >= rr2 .or. rr2 < eps_fp .or. kcs(nm)==3 .or. kcs(num)==3) then
                            r00 = r0(nm, k, l)
                         else
                            r00 = r0(nm , k, l)                                     &
                                & + (1.0_fp - cfl)*(r0(nm , k, l) - r0(ndm, k, l))  &
                                &                 *(r0(num, k, l) - r0(nm , k, l))  &
                                &                 /(r0(num, k, l) - r0(ndm, k, l))
                         endif
                         flux = qyk(nm, k) * r00
                         if (kcs(nm) == 1) then
                            fluxv(nm, k, l) = fluxv(nm, k, l) + flux
                         endif
                      enddo
                   elseif (qyk(nm, k) < 0.0_fp) then
                      do l = 1, lstsci
                         rr1 = abs(r0(nuum, k, l) - 2.0_fp*r0(num, k, l) + r0(nm, k, l))
                         rr2 = abs(r0(nuum, k, l) - r0(nm, k, l))
                         if (kfsz0(nm, k)*kfsz0(nuum, k) == 0 .or. rr1 >= rr2 .or. rr2 < eps_fp .or. kcs(nm)==3 .or. kcs(num)==3) then
                            r00 = r0(num, k, l)
                         else
                            r00 = r0(num, k, l)                                      &
                                & + (1.0_fp + cfl)*(r0(nm , k, l) - r0(num , k, l))  &
                                &                 *(r0(num, k, l) - r0(nuum, k, l))  &
                                &                 /(r0(nm , k, l) - r0(nuum, k, l))
                         endif
                         flux = qyk(nm, k) * r00
                         if (kcs(nm) == 1) then
                            fluxv(nm, k, l) = fluxv(nm, k, l) + flux
                         endif
                      enddo
                   endif
                endif
             enddo
          endif
          !
          ! DIFFUSIVE TRANSPORT IN Y-DIRECTION
          !
          if (kfs(nm)*kfs(num) == 1) then
             do k = kfvmin(nm), kmax
                if (kfvz0(nm, k)==1 .and. kfsz0(nm, k)*kfsz0(num, k)==1) then
                   difl = dicuv(nm, k)
                   difr = dicuv(num, k)
                   kenv = max(0, abs(2 - kcs(nm)))
                   do l = 1, lstsci
                      cl              = r0(nm, k, l)
                      cr              = r0(num, k, l)
                      flux            = 0.5_fp*(cr - cl)*(difl + difr)/(sigdif(l)*guv(nm))
                      fluxv(nm, k, l) = fluxv(nm, k, l) - areav(nm, k)*flux*kenv
                   enddo
                endif
             enddo
          endif
       enddo

    elseif(trasol == 'iupw         ') then

       do nm = 1, nmmax
          !
          ! ADVECTIVE TRANSPORT IN X-DIRECTION
          !
          nmu  = nm + icx
          num  = nm + icy
          if (kfu(nm) /= 0 .and. kcs(nm)*kcs(nmu) /= 0) then
             kmin   = max(kfumin(nm), 1)
             kkmin  = min(kfsmx0(nm), kfsmax(nm))
             kkminu = min(kfsmx0(nmu), kfsmax(nmu))
             do k = kmin, min(kkmin, kkminu)
                do l = 1, lstsci
                   if (qxk(nm,k) > 0.0_fp) then
                      fluxu(nm,k,l) = fluxu(nm,k,l) + qxk(nm, k) * r1(nm , k, l)
                   elseif (qxk(nm,k) < 0.0_fp) then
                      fluxu(nm,k,l) = fluxu(nm,k,l) + qxk(nm, k) * r1(nmu, k, l)
                   endif
                enddo
             enddo
          endif
          !
          ! Adjacent cell water fall correction (explicit)
          !
          if (kfumx0(nm) > min(kkmin,kkminu)) then
             do k = min(kkmin,kkminu)+1, kfumx0(nm)
                if (qxk(nm, k) > 0.0_fp) then
                   do l = 1, lstsci
                      flux = qxk(nm, k) * r0(nm , k, l)
                      if (kcs(nm) == 1) then
                         fluxu(nm, k, l) = fluxu(nm, k, l) + flux
                      endif
                   enddo
                elseif (qxk(nm, k) < 0.0_fp) then
                   do l = 1, lstsci
                      flux = qxk(nm, k) * r0(num, k, l)
                      if (kcs(nm) == 1) then
                         fluxu(nm, kkmin, l) = fluxu(nm, kkmin, l) + flux
                      endif
                   enddo
                endif
             enddo
          endif
          !
          ! ADVECTIVE TRANSPORT IN Y-DIRECTION
          !
          if (kfv(nm) /= 0  .and. kcs(nm)*kcs(num) /= 0) then
             kmin   = max(kfvmin(nm), 1)
             kkmin  = min(kfsmx0(nm), kfsmax(nm))
             kkminv = min(kfsmx0(num), kfsmax(num))
             do k = kmin, min(kkmin,kkminv)
                do l = 1, lstsci
                   if (qyk(nm,k) > 0.0_fp) then
                      fluxv(nm,k,l) = fluxv(nm,k,l) + qyk(nm, k) * r1(nm , k, l)
                   elseif (qyk(nm,k) < 0.0_fp) then
                      fluxv(nm,k,l) = fluxv(nm,k,l) + qyk(nm, k) * r1(num, k, l)
                   endif
                enddo
             enddo
          endif
          !
          ! Adjacent cell water fall correction (explicit)
          !
          if (kfvmx0(nm) > min(kkmin,kkminv)) then
             do k = min(kkmin,kkminv)+1, kfvmx0(nm)
                if (qyk(nm, k) > 0.0_fp) then
                   do l = 1, lstsci
                      flux = qyk(nm, k) * r0(nm , k, l)
                      if (kcs(nm) == 1) then
                         fluxv(nm, k, l) = fluxv(nm, k, l) + flux
                      endif
                   enddo
                elseif (qyk(nm, k) < 0.0_fp) then
                   do l = 1, lstsci
                      flux = qyk(nm, k) * r0(num, k, l)
                      if (kcs(nm) == 1) then
                         fluxv(nm, kkmin, l) = fluxv(nm, kkmin, l) + flux
                      endif
                   enddo
                endif
             enddo
          endif
          !
          ! DIFFUSIVE TRANSPORT IN X-DIRECTION
          !
          if (kfu(nm) == 1) then
             mink_old = min(kfsmx0(nm),kfsmx0(nmu))
             mink_new = min(kfsmax(nm),kfsmax(nmu))
             mink     = min(mink_old  ,mink_new)
             do k = kfumin(nm), mink
                if (kfuz0(nm, k)==1 .and. kfsz0(nm, k)*kfsz0(nmu, k)==1) then
                   difl = dicuv(nm, k)
                   difr = dicuv(nmu, k)
                   flux = 0.5_fp * (difl + difr) / (0.7_fp*gvu(nm))
                   kenu = max(0, abs(2 - kcs(nm)))
                   do l = 1, lstsci
                      fluxu(nm,k,l) = fluxu(nm,k,l) + areau(nm, k)*flux*real(kenu,fp)
                   enddo
                endif
             enddo
          endif
          !
          ! DIFFUSIVE TRANSPORT IN Y-DIRECTION
          !
          if (kfv(nm) == 1) then
             mink_old = min(kfsmx0(nm),kfsmx0(num))
             mink_new = min(kfsmax(nm),kfsmax(num))
             mink     = min(mink_old  ,mink_new)
             do k = kfvmin(nm), mink
                if (kfvz0(nm, k)==1 .and. kfsz0(nm, k)*kfsz0(num, k)==1) then
                   difl = dicuv(nm, k)
                   difr = dicuv(num, k)
                   flux = 0.5_fp * (difl + difr) / (0.7_fp*guv(nm))
                   kenv = max(0, abs(2 - kcs(nm)))
                   do l = 1, lstsci
                      fluxv(nm,k,l) = fluxv(nm,k,l) + areav(nm, k)*flux*real(kenv,fp)
                   enddo
                endif
             enddo
          endif
       enddo
    endif
    !
    ! Cumulative flux
    !
    if (flwoutput%cumdifuflux) then
       if (icx == 1) then
         fluxuc = fluxuc + fluxu * timest
         fluxvc = fluxvc + fluxv * timest
       else
         fluxuc = fluxuc + fluxv * timest
         fluxvc = fluxvc + fluxu * timest
       endif
    endif
    !
    ! swap fluxu and fluxv for sediment transport in z layer morphology.
    ! Originally, switch was a full sized array and the code was without do-loops:
    !     switch = fluxu
    !     fluxu  = fluxv
    !     fluxv  = fluxu
    ! These statements may cause problems when the arrays are very big; values are put
    ! on the (restricted) stack
    ! Therefore, the primitive triple do loop is used, where switch can be a scalar
    !
    if (icx == 1) then
       do nm =  gdp%d%nmlb,gdp%d%nmub
           do k = 1, kmax
               do l =  1,lstsci
                 switch          = fluxu (nm,k,l)
                 fluxu (nm,k,l)  = fluxv (nm,k,l)
                 fluxv (nm,k,l)  = switch
               enddo
           enddo
       enddo
    endif
end subroutine z_difuflux

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
module test
  contains
subroutine printerror(host, errorcode, privatedata) BIND(C)
  use ISO_C_BINDING
  use OPI
  type(c_ptr), value :: host
  integer(c_int), value :: errorcode
  type(c_ptr), value :: privatedata
  write(*,*) 'OPI Error: ', OPI_ErrorMessage(errorcode)

end subroutine
end module
program opi_host_fortran_test
  use OPI
  use OPI_Types
  use test

  type(c_ptr) :: host
  type(c_ptr) :: data
  type(c_ptr) :: propagator
  integer(C_INT) :: status
  type(OPI_Orbit),pointer :: orbit (:)
  integer(C_INT) :: testint = 10
  type(c_ptr) :: mode
  ! create host
  host = OPI_createHost()

  call OPI_setErrorCallback(host, C_FUNLOC(printerror), C_NULL_PTR)

  ! load plugins
  status = OPI_loadPlugins(host, "plugins")

  ! create data object AFTER! loading all support plugins
  data = OPI_createData(host, 200)

  write(*,*) 'Registered Propagator:'
  do i = 1, OPI_getPropagatorCount(host), 1
    write(*,*) '# ', i, ': ', OPI_Module_getName(OPI_getPropagator(host, i - 1))
  end do

  ! find a specific propagator
  propagator = OPI_getPropagator(host, "CPP Example Propagator")

  ! a null-pointer (0) is returned if no propagator is found
  if(.NOT. C_ASSOCIATED(propagator)) then
    write(*,*) 'Propagator not found'
  else
    write(*,*) 'Propagator Properties:'
    do i= 1, OPI_Module_getPropertyCount(propagator), 1
      write(*,*) '# ', i, ': ', OPI_Module_getPropertyName(propagator, i -1), " value: ", &
      OPI_Module_getPropertyString(propagator, i -1)
    end do
    ! run propagation of our test propagator
    status = OPI_Propagator_propagate(propagator, data, 0.d0, 0.d0, mode, c_null_ptr)

    ! retrieve orbital data values
    orbit => OPI_Population_getOrbit(data)
    ! print inclination values
    do i = 1, 200, 1
      write(*,*) orbit(i)%inclination
    end do
  end if

  ! free data
  status = OPI_destroyData(data)

  ! drop host
  status = OPI_destroyHost(host)
end program opi_host_fortran_test

max_stars_repo_name: Space-Systems/OPI

max_issues_repo_name: Space-Systems/OPI

max_forks_repo_name: Space-Systems/OPI


text:
c
c
c     ############################################################
c     ##  COPYRIGHT (C) 1995 by Yong Kong & Jay William Ponder  ##
c     ##                  All Rights Reserved                   ##
c     ############################################################
c
c     #################################################################
c     ##                                                             ##
c     ##  subroutine rotpole  --  rotate multipoles to global frame  ##
c     ##                                                             ##
c     #################################################################
c
c
c     "rotpole" constructs the set of atomic multipoles in the global
c     frame by applying the correct rotation matrix for each site
c
c
      subroutine rotpole
      implicit none
      include 'sizes.i'
      include 'mpole.i'
      integer i
      real*8 a(3,3)
c
c
c     rotate the atomic multipoles at each site in turn
c
      do i = 1, npole
         call rotmat (i,a)
         call rotsite (i,a)
      end do
      return
      end
c
c
c     ################################################################
c     ##                                                            ##
c     ##  subroutine rotmat  --  find global frame rotation matrix  ##
c     ##                                                            ##
c     ################################################################
c
c
c     "rotmat" finds the rotation matrix that converts from the local
c     coordinate system to the global frame at a multipole site
c
c
      subroutine rotmat (i,a)
      implicit none
      include 'sizes.i'
      include 'atoms.i'
      include 'mpole.i'
      integer i,ii
      integer ix,iy,iz
      real*8 r,dot
      real*8 random
      real*8 xi,yi,zi
      real*8 dx,dy,dz
      real*8 dx1,dy1,dz1
      real*8 dx2,dy2,dz2
      real*8 dx3,dy3,dz3
      real*8 a(3,3)
c
c
c     get coordinates and frame definition for the multipole site
c
      ii = ipole(i)
      xi = x(ii)
      yi = y(ii)
      zi = z(ii)
      ix = xaxis(i)
      iy = yaxis(i)
      iz = zaxis(i)
c
c     use the identity matrix as the default rotation matrix
c
      a(1,1) = 1.0d0
      a(2,1) = 0.0d0
      a(3,1) = 0.0d0
      a(1,3) = 0.0d0
      a(2,3) = 0.0d0
      a(3,3) = 1.0d0
c
c     Z-Only method rotation matrix elements for z-axis only
c
      if (polaxe(i) .eq. 'Z-Only') then
         dx = x(iz) - xi
         dy = y(iz) - yi
         dz = z(iz) - zi
         r = sqrt(dx*dx + dy*dy + dz*dz)
         a(1,3) = dx / r
         a(2,3) = dy / r
         a(3,3) = dz / r
         dx = random ()
         dy = random ()
         dz = random ()
         dot = dx*a(1,3) + dy*a(2,3) + dz*a(3,3)
         dx = dx - dot*a(1,3)
         dy = dy - dot*a(2,3)
         dz = dz - dot*a(3,3)
         r = sqrt(dx*dx + dy*dy + dz*dz)
         a(1,1) = dx / r
         a(2,1) = dy / r
         a(3,1) = dz / r
c
c     Z-then-X method rotation matrix elements for z- and x-axes
c
      else if (polaxe(i) .eq. 'Z-then-X') then
         dx = x(iz) - xi
         dy = y(iz) - yi
         dz = z(iz) - zi
         r = sqrt(dx*dx + dy*dy + dz*dz)
         a(1,3) = dx / r
         a(2,3) = dy / r
         a(3,3) = dz / r
         dx = x(ix) - xi
         dy = y(ix) - yi
         dz = z(ix) - zi
         dot = dx*a(1,3) + dy*a(2,3) + dz*a(3,3)
         dx = dx - dot*a(1,3)
         dy = dy - dot*a(2,3)
         dz = dz - dot*a(3,3)
         r = sqrt(dx*dx + dy*dy + dz*dz)
         a(1,1) = dx / r
         a(2,1) = dy / r
         a(3,1) = dz / r
c
c     Bisector method rotation matrix elements for z- and x-axes
c
      else if (polaxe(i) .eq. 'Bisector') then
         dx = x(iz) - xi
         dy = y(iz) - yi
         dz = z(iz) - zi
         r = sqrt(dx*dx + dy*dy + dz*dz)
         dx1 = dx / r
         dy1 = dy / r
         dz1 = dz / r
         dx = x(ix) - xi
         dy = y(ix) - yi
         dz = z(ix) - zi
         r = sqrt(dx*dx + dy*dy + dz*dz)
         dx2 = dx / r
         dy2 = dy / r
         dz2 = dz / r
         dx = dx1 + dx2
         dy = dy1 + dy2
         dz = dz1 + dz2
         r = sqrt(dx*dx + dy*dy + dz*dz)
         a(1,3) = dx / r
         a(2,3) = dy / r
         a(3,3) = dz / r
         dot = dx2*a(1,3) + dy2*a(2,3) + dz2*a(3,3)
         dx = dx2 - dot*a(1,3)
         dy = dy2 - dot*a(2,3)
         dz = dz2 - dot*a(3,3)
         r = sqrt(dx*dx + dy*dy + dz*dz)
         a(1,1) = dx / r
         a(2,1) = dy / r
         a(3,1) = dz / r
c
c     Z-Bisect method rotation matrix elements for z- and x-axes
c
      else if (polaxe(i) .eq. 'Z-Bisect') then
         dx = x(iz) - xi
         dy = y(iz) - yi
         dz = z(iz) - zi
         r = sqrt(dx*dx + dy*dy + dz*dz)
         a(1,3) = dx / r
         a(2,3) = dy / r
         a(3,3) = dz / r
         dx = x(ix) - xi
         dy = y(ix) - yi
         dz = z(ix) - zi
         r = sqrt(dx*dx + dy*dy + dz*dz)
         dx1 = dx / r
         dy1 = dy / r
         dz1 = dz / r
         dx = x(iy) - xi
         dy = y(iy) - yi
         dz = z(iy) - zi
         r = sqrt(dx*dx + dy*dy + dz*dz)
         dx2 = dx / r
         dy2 = dy / r
         dz2 = dz / r
         dx = dx1 + dx2
         dy = dy1 + dy2
         dz = dz1 + dz2
         r = sqrt(dx*dx + dy*dy + dz*dz)
         dx = dx / r
         dy = dy / r
         dz = dz / r
         dot = dx*a(1,3) + dy*a(2,3) + dz*a(3,3)
         dx = dx - dot*a(1,3)
         dy = dy - dot*a(2,3)
         dz = dz - dot*a(3,3)
         r = sqrt(dx*dx + dy*dy + dz*dz)
         a(1,1) = dx / r
         a(2,1) = dy / r
         a(3,1) = dz / r
c
c     3-Fold method rotation matrix elements for z- and x-axes
c
      else if (polaxe(i) .eq. '3-Fold') then
         dx = x(iz) - xi
         dy = y(iz) - yi
         dz = z(iz) - zi
         r = sqrt(dx*dx + dy*dy + dz*dz)
         dx1 = dx / r
         dy1 = dy / r
         dz1 = dz / r
         dx = x(ix) - xi
         dy = y(ix) - yi
         dz = z(ix) - zi
         r = sqrt(dx*dx + dy*dy + dz*dz)
         dx2 = dx / r
         dy2 = dy / r
         dz2 = dz / r
         dx = x(iy) - xi
         dy = y(iy) - yi
         dz = z(iy) - zi
         r = sqrt(dx*dx + dy*dy + dz*dz)
         dx3 = dx / r
         dy3 = dy / r
         dz3 = dz / r
         dx = dx1 + dx2 + dx3
         dy = dy1 + dy2 + dy3
         dz = dz1 + dz2 + dz3
         r = sqrt(dx*dx + dy*dy + dz*dz)
         a(1,3) = dx / r
         a(2,3) = dy / r
         a(3,3) = dz / r
         dot = dx2*a(1,3) + dy2*a(2,3) + dz2*a(3,3)
         dx = dx2 - dot*a(1,3)
         dy = dy2 - dot*a(2,3)
         dz = dz2 - dot*a(3,3)
         r = sqrt(dx*dx + dy*dy + dz*dz)
         a(1,1) = dx / r
         a(2,1) = dy / r
         a(3,1) = dz / r
      end if
c
c     finally, find rotation matrix elements for the y-axis
c
      a(1,2) = a(3,1)*a(2,3) - a(2,1)*a(3,3)
      a(2,2) = a(1,1)*a(3,3) - a(3,1)*a(1,3)
      a(3,2) = a(2,1)*a(1,3) - a(1,1)*a(2,3)
      return
      end
c
c
c     ################################################################
c     ##                                                            ##
c     ##  subroutine rotsite  --  rotate multipoles at single site  ##
c     ##                                                            ##
c     ################################################################
c
c
c     "rotsite" computes the atomic multipoles at a specified site
c     in the global coordinate frame by applying a rotation matrix
c
c
      subroutine rotsite (isite,a)
      implicit none
      include 'sizes.i'
      include 'atoms.i'
      include 'mpole.i'
      integer i,j,k,m
      integer isite
      real*8 a(3,3)
      real*8 m2(3,3)
      real*8 r2(3,3)
c
c
c     monopoles have the same value in any coordinate frame
c
      rpole(1,isite) = pole(1,isite)
c
c     rotate the dipoles to the global coordinate frame
c
      do i = 2, 4
         rpole(i,isite) = 0.0d0
         do j = 2, 4
            rpole(i,isite) = rpole(i,isite) + pole(j,isite)*a(i-1,j-1)
         end do
      end do
c
c     rotate the quadrupoles to the global coordinate frame
c
      k = 5
      do i = 1, 3
         do j = 1, 3
            m2(i,j) = pole(k,isite)
            r2(i,j) = 0.0d0
            k = k + 1
         end do
      end do
      do i = 1, 3
         do j = 1, 3
            if (j .lt. i) then
               r2(i,j) = r2(j,i)
            else
               do k = 1, 3
                  do m = 1, 3
                     r2(i,j) = r2(i,j) + a(i,k)*a(j,m)*m2(k,m)
                  end do
               end do
            end if
         end do
      end do
      k = 5
      do i = 1, 3
         do j = 1, 3
            rpole(k,isite) = r2(i,j)
            k = k + 1
         end do
      end do
      return
      end

ext: f

max_stars_repo_name: htna/HCsbLib

max_issues_repo_name: htna/HCsbLib

max_forks_repo_name: htna/HCsbLib


max_stars_repo_name: manoelfranca/cilppp

max_issues_repo_name: manoelfranca/cilppp

max_forks_repo_name: manoelfranca/cilppp


text: #redirect Uncle Vitos Slice of N.Y.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


!> \file mo_spatialsimilarity.f90

!> \brief Routines for bias insensitive comparison of spatial patterns.

!> \details These routines are based on the idea that spatial similarity can be assessed by comparing
!>          the magnitude of neighboring pixels (e.g. is the neighboring pixel larger or smaller).

!> \author Matthias Zink, Juliane Mai
!> \date Mar 2013

MODULE mo_spatialsimilarity

  ! This module contains routines for the masked calculation of
  ! spatial similarity measures like NDV

  ! Written Nov 2012, Matthias Zink and Juliane Mai

  ! License
  ! -------
  ! This file is part of the JAMS Fortran package, distributed under the MIT License.
  !
  ! Copyright (c) 2012-2015 Matthias Zink, Juliane Mai
  !
  ! Permission is hereby granted, free of charge, to any person obtaining a copy
  ! of this software and associated documentation files (the "Software"), to deal
  ! in the Software without restriction, including without limitation the rights
  ! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
  ! copies of the Software, and to permit persons to whom the Software is
  ! furnished to do so, subject to the following conditions:
  !
  ! The above copyright notice and this permission notice shall be included in all
  ! copies or substantial portions of the Software.
  !
  ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
  ! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  ! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
  ! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  ! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
  ! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  ! SOFTWARE.

  USE mo_kind,   ONLY: i4, sp, dp

  IMPLICIT NONE

  PUBLIC :: NNDV       ! number of neighboring dominating values
  PUBLIC :: PD         ! patter dissimilarity measure

  ! ------------------------------------------------------------------

  !     NAME
  !         NNDV - number of neighboring dominating values

  !     PURPOSE
  !>         \brief    Calculates the number of neighboring dominating values, a measure for spatial dissimilarity.
  !>         \details
  !>             NNDV             = 1 - sum(abs(dominating_neighbors(mat1) - dominating_neighbors(mat2))) / count(mask)
  !>             dominating_neighbors(mat1) = comparison if pixel is larger than its neighbouring values
  !>
  !>            An array element value is compared with its 8 neighbouring cells to check if these cells are larger
  !>            than the array element value. The result is a 3x3 matrix in which larger cells are indicated by a
  !>            true value. For comparison this is done with both input arrays. The resulting array is the sum of
  !>            substraction of the 3x3 matrices for each of the both arrays. The resulting matrix is afterwards
  !>            normalized to its available neighbors.
  !>            Furthermore an average over the entire field is calculated. The valid interval of
  !>            the values for NNDV is [0..1]. In which 1 indicates full agreement and 0 full dismatching.
  !>

  !>            EXAMPLE:\n
  !>            mat1 =  | 12 17  1 | , mat2 = |  7  9 12 |
  !>                    |  4 10 11 |          | 12 11 11 |
  !>                    | 15  2 20 |          |  5 13  7 |
  !>            booleans determined for every grid cell following fortran array scrolling
  !>            i.e. (/col1_row1, col1_row2, col1_row3, col2_row1, .. ,col3_row3/),(/3,3/)
  !>
  !>            comp1 = | FFF FFF FTF, FFF FFF FFF, FTT FFT FFF |
  !>                    | FFF TFT TTF, TFT TFF FTT, TFF FFT FFF |
  !>                    | FFF FFF FFF, TTF TFF TTF, FFF FFF FFF |
  !>
  !>            comp2 = | FFF FFT FTT, FFT FFT FTT, FFF FFF FFF |
  !>                    | FFF FFF FFT, FTF FFT TFF, FFT TFF FFF |
  !>                    | FFF TFF TTF, FFF FFF FFF, TTF TFF FFF |
  !>
  !>           NNDVMatrix =
  !>           abs( count(comp1) - count(comp2) ) = | 1-3, 0-4, 3-0 | = | 2, 4, 3 |
  !>                                                | 4-1, 5-3, 2-2 |   | 3, 2, 0 |
  !>                                                | 0-3, 5-0, 0-3 |   | 3, 5, 3 |
  !>
  !>                                    DISSIMILAR / VALID NEIGH CELLS
  !>           NNDVMatrix / VALID NEIGH CELLS = | 2, 4, 3 | / | 3, 5, 3 |
  !>                                            | 3, 2, 0 |   | 5, 8, 5 |
  !>                                            | 3, 5, 3 |   | 3, 5, 3 |
  !>
  !>                                          = | 0.66, 0.80, 1.00 |
  !>                                            | 0.60, 0.25, 0.00 |
  !>                                            | 1.0,  1.00, 1.00 |
  !>
  !>           NNDV = 1 - sum(NNDVMatrix) / count(mask) = 1 - (6.31666666 / 9) = 0.2981
  !>

  !>
  !>           If an optinal mask is given, the calculations are over those locations that correspond to true values in the mask.
  !>           mat1 and mat2 can be single or double precision. The result will have the same numerical precision.

  !     CALLING SEQUENCE
  !         out = NNDV(mat1, mat2, mask=mask, valid=valid)

  !     INDENT(IN)
  !>        \param[in] "real(sp/dp), dimension(:,:) :: mat1" 2D-array with input numbers
  !>        \param[in] "real(sp/dp), dimension(:,:) :: mat2" 2D-array with input numbers

  !     INDENT(INOUT)
  !         None

  !     INDENT(OUT)
  !        None

  !     INDENT(IN), OPTIONAL
  !>        \param[in] "logical,dimension(:,:),optinal :: mask" 2D-array of logical values with size(mat1/mat2).
  !>           If present, only those locations in mat1/mat2 having true values in mask are evaluated.

  !     INDENT(INOUT), OPTIONAL
  !         None

  !     INDENT(OUT), OPTIONAL
  !>        \param[out] "logical              :: valid"   indicates if the function could determine a valid value
  !>                                                      result can be unvalid if entire mask is .false. for ex.
  !>                                                      in this case PatternDissim is set to 0 (worst case)

  !     RETURN
  !>        \return real(sp/dp) :: NNDV — Number of neighboring dominating values

  !     RESTRICTIONS
  !         Input values must be floating points.

  !     EXAMPLE
  !         mat1 = reshape(/ 1., 2, 3., -999., 5., 6. /, (/3,3/))
  !         mat2 = reshape(/ 1., 2, 3., -999., 5., 6. /, (/3,3/))
  !         out  = NNDV(mat1, mat2, mask=(mat1 >= 0. .and. mat2 >= 0.))
  !         -> see also example in test directory

  !     LITERATURE
  !>         \note routine based on algorithm by Luis Samaniego 2009

  !     HISTORY
  !>         \author Matthias Zink
  !>         \date   Nov 2012
  !          update  May 2015 created documentation
  INTERFACE NNDV
     MODULE PROCEDURE NNDV_sp, NNDV_dp
  END INTERFACE NNDV

  ! ------------------------------------------------------------------

  !     NAME
  !         PD

  !     PURPOSE
  !>         \brief Calculates pattern dissimilarity (PD) measure
  !>         \details
  !>             PD             = 1 - sum(dissimilarity(mat1, mat2)) / count(mask)
  !>             dissimilarity(mat1, mat2) = comparison if pixel is larger than its neighbouring values
  !>
  !>            An array element value is compared with its 8 neighbouring cells to check if these cells are larger
  !>            than the array element value. The result is a 3x3 matrix in which larger cells are indicated by a
  !>            true value. For comparison this is done with both input arrays. The resulting array is the sum of
  !>            xor values of the 3x3 matrices for each of the both arrays.  This means only neighbourhood comparisons
  !>            which are different in the 2 matrices are counted. This resulting matrix is afterwards normalized to its
  !>            available neighbors. Furthermore an average over the entire field is calculated. The valid interval of
  !>            the values for PD is [0..1]. In which 1 indicates full agreement and 0 full dismatching.
  !>
  !>

  !>            EXAMPLE:\n
  !>            mat1 =  | 12 17  1 | , mat2 = |  7  9 12 |
  !>                    |  4 10 11 |          | 12 11 11 |
  !>                    | 15  2 20 |          |  5 13  7 |
  !>
  !>            booleans determined for every grid cell following fortran array scrolling
  !>            i.e. (/col1_row1, col1_row2, col1_row3, col2_row1, .. ,col3_row3/),(/3,3/)
  !>
  !>            comp1 = | FFF FFF FTF, FFF FFF FFF, FTT FFT FFF |
  !>                    | FFF TFT TTF, TFT TFF FTT, TFF FFT FFF |
  !>                    | FFF FFF FFF, TTF TFF TTF, FFF FFF FFF |
  !>
  !>            comp2 = | FFF FFT FTT, FFT FFT FTT, FFF FFF FFF |
  !>                    | FFF FFF FFT, FTF FFT TFF, FFT TFF FFF |
  !>                    | FFF TFF TTF, FFF FFF FFF, TTF TFF FFF |
  !>
  !>
  !>            xor=neq = | FFF FFT FFT, FFT FFT FTT, FTT FFT FFF |
  !>                      | FFF TFT TTT, TTT TFT TTT, TFT TFT FFF |
  !>                      | FFF TFF TTF, TTF TFF TTF, TTF TFF FFF |
  !>
  !>                        DISSIMILAR / VALID NEIGH CELLS
  !>            PDMatrix = | 2, 4, 3 | / | 3, 5, 3 | = | 0.66, 0.80, 1.00 |
  !>                       | 5, 8, 4 |   | 5, 8, 5 |   | 1.00, 1.00, 0.80 |
  !>                       | 3, 5, 3 |   | 3, 5, 3 |   | 1.00, 1.00, 1.00 |
  !>
  !>           PD = 1 - sum(PDMatrix) / count(mask) = 1 - (8.2666666 / 9) = 0.08148
  !>

  !>
  !>           If an optinal mask is given, the calculations are over those locations that correspond to true values in the mask.
  !>           mat1 and mat2 can be single or double precision. The result will have the same numerical precision.

  !     CALLING SEQUENCE
  !         out = PD(mat1, mat2, mask=mask, valid=valid)

  !     INDENT(IN)
  !>        \param[in] "real(sp/dp), dimension(:,:) :: mat1" 2D-array with input numbers
  !>        \param[in] "real(sp/dp), dimension(:,:) :: mat2" 2D-array with input numbers

  !     INDENT(INOUT)
  !         None

  !     INDENT(OUT)
  !         None

  !     INDENT(IN), OPTIONAL
  !>        \param[in] "logical,dimension(:,:),optinal :: mask" 2D-array of logical values with size(mat1/mat2)
  !>           If present, only those locations in mat1/mat2 having true values in mask are evaluated.

  !     INDENT(INOUT), OPTIONAL
  !         None

  !     INDENT(OUT), OPTIONAL
  !>        \param[out] "logical,optinal :: valid"   indicates if the function could determine a valid value
  !>                                                      result can be unvalid if entire mask is .false. for ex.
  !>                                                      in this case PD is set to 0 (worst case)

  !     RETURN
  !>        \return real(sp/dp) :: PD — pattern dissimilarity measure

  !     RESTRICTIONS
  !         Input values must be floating points.

  !     EXAMPLE
  !         mat1 = reshape(/ 1., 2, 3., -999., 5., 6. /, (/3,3/))
  !         mat2 = reshape(/ 1., 2, 3., -999., 5., 6. /, (/3,3/))
  !         out  = PD(mat1, mat2, mask=(mat1 >= 0. .and. mat2 >= 0.))
  !         -> see also example in test directory

  !     LITERATURE
  !          None

  !     HISTORY
  !>         \author Matthias Zink and Juliane Mai
  !>         \date   Jan 2013
  INTERFACE PD
     MODULE PROCEDURE PD_sp, PD_dp
  END INTERFACE PD

  ! ------------------------------------------------------------------

CONTAINS

  FUNCTION NNDV_sp(mat1, mat2, mask, valid)

    IMPLICIT NONE

    REAL(sp),    DIMENSION(:,:),                                            INTENT(IN)  :: mat1, mat2
    LOGICAL,     DIMENSION(:,:),                                  OPTIONAL, INTENT(IN)  :: mask
    LOGICAL,                                                      OPTIONAL, INTENT(OUT) :: valid
    REAL(sp)                                                                            :: NNDV_sp

    INTEGER(i4)                                                                         :: iCo, iRo
    INTEGER(i4)                                                                         :: noValidPixels
    INTEGER(i4), DIMENSION(size(shape(mat1)) )                                          :: shapemask
    INTEGER(i4), DIMENSION(size(mat1, dim=1), size(mat1, dim=2))                        :: validcount
    REAL(sp),    DIMENSION(size(mat1, dim=1)+2_i4, size(mat1, dim=2)+2_i4)              :: bufferedMat1, bufferedMat2
    REAL(sp)   , DIMENSION(size(mat1, dim=1), size(mat1, dim=2))                        :: NNDVMatrix
    LOGICAL,     DIMENSION(size(mat1, dim=1)+2_i4, size(mat1, dim=2)+2_i4)              :: maske
    LOGICAL,     DIMENSION(size(mat1, dim=1), size(mat1, dim=2))                        :: validmaske

    ! check if input has all the same dimensions
    if (present(mask)) then
       shapemask = shape(mask)
    else
       shapemask = shape(mat1)
    end if
    !
    if (any(shape(mat1) .NE. shape(mat2)))  &
         stop 'NNDV_sp: shapes of input matrix 1 and input matrix 2 are not matching'
    if (any(shape(mat1) .NE. shapemask))    &
         stop 'NNDV_sp: shapes of input matrices and mask are not matching'
    !
    ! additional 2 rows and 2 cols added for checking the border cells without crashing the search agorithm
    ! so the search windows can always cover 9 cells even if it is checking a border cell
    ! buffer rows are initialized as false values within mask, i.e. they are not considered for the calculation
    ! of the criteria
    !
    ! initialize mask with default=.false.
    maske = .false.
    if (present(mask)) then
       maske(2:(size(maske,dim=1)-1_i4), 2:(size(maske,dim=2)-1_i4)) = mask
    else
       maske(2:(size(maske,dim=1)-1_i4), 2:(size(maske,dim=2)-1_i4)) = .true.
    end if
    !
    ! initialize bufferedMat1 & bufferedMat2
    bufferedMat1 = 0.0_sp
    bufferedMat2 = 0.0_sp
    bufferedMat1(2:(size(maske,dim=1)-1), 2:(size(maske,dim=2)-1)) = mat1
    bufferedMat2(2:(size(maske,dim=1)-1), 2:(size(maske,dim=2)-1)) = mat2
    !
    ! initialize NNDV
    NNDVMatrix = 0.0_sp
    !
    NNDV_sp = 0.0_sp
    do iCo = 2_i4, size(bufferedMat1, dim=2) - 1
       do iRo = 2_i4, size(bufferedMat1, dim=1) - 1
          if (.NOT. maske(iRo,iCo)) cycle
          NNDVMatrix(iRo-1_i4,iCo-1_i4) =   &
               real(    &
                  abs(  &
                     count((bufferedMat1(iRo-1:iRo+1 , iCo-1:iCo+1) - bufferedMat1(iRo,iCo) > epsilon(0.0_sp)) .AND. &
                           (maske(iRo-1:iRo+1 , iCo-1:iCo+1))) -                                                     &
                     count((bufferedMat2(iRo-1:iRo+1 , iCo-1:iCo+1) - bufferedMat2(iRo,iCo) > epsilon(0.0_sp)) .AND. &
                           (maske(iRo-1:iRo+1 , iCo-1:iCo+1)))                                                       &
                  ),    &
               sp)
          ! count - 1 to exclude referendce cell (iRo, iCo)
          validcount(iRo-1_i4,iCo-1_i4) = count(maske(iRo-1_i4:iRo+1_i4 , iCo-1_i4:iCo+1_i4)) - 1_i4
        end do
    end do
    !
    ! normalize every pixel to number of valid neighbouring cells (defined by maske) [0..8] --> [0..1]
    validmaske = (maske(2:size(maske, dim=1) - 1_i4, 2:size(maske, dim=2) - 1_i4) .and. (validcount > 0_i4))
    noValidPixels = count(validmaske)
    if (noValidPixels .GT. 0_i4) then
       NNDVMatrix = merge(NNDVMatrix/real(validcount, sp),NNDVMatrix, validmaske)
       ! average over all pixels
       NNDV_sp = 1.0_sp - sum(NNDVMatrix, mask=validmaske) / noValidPixels
       if (present(valid)) valid = .TRUE.
    ! case if maske is full of .false. or no valid neighbouring cells available for every pixel (validcount(:,:)=0)
    else
       NNDV_sp = 0.0_sp
       if (present(valid)) valid = .FALSE.
    end if
    !
  END FUNCTION NNDV_sp

 FUNCTION NNDV_dp(mat1, mat2, mask, valid)

    IMPLICIT NONE

    REAL(dp),    DIMENSION(:,:),                                            INTENT(IN)  :: mat1, mat2
    LOGICAL,     DIMENSION(:,:),                                  OPTIONAL, INTENT(IN)  :: mask
    LOGICAL,                                                      OPTIONAL, INTENT(OUT) :: valid
    REAL(dp)                                                                            :: NNDV_dp

    INTEGER(i4)                                                                         :: iCo, iRo
    INTEGER(i4)                                                                         :: noValidPixels
    INTEGER(i4), DIMENSION(size(shape(mat1)) )                                          :: shapemask
    INTEGER(i4), DIMENSION(size(mat1, dim=1), size(mat1, dim=2))                        :: validcount
    REAL(dp),    DIMENSION(size(mat1, dim=1)+2_i4, size(mat1, dim=2)+2_i4)              :: bufferedMat1, bufferedMat2
    REAL(dp)   , DIMENSION(size(mat1, dim=1), size(mat1, dim=2))                        :: NNDVMatrix
    LOGICAL,     DIMENSION(size(mat1, dim=1)+2_i4, size(mat1, dim=2)+2_i4)              :: maske
    LOGICAL,     DIMENSION(size(mat1, dim=1), size(mat1, dim=2))                        :: validmaske

    ! check if input has all the same dimensions
    if (present(mask)) then
       shapemask = shape(mask)
    else
       shapemask = shape(mat1)
    end if
    !
    if (any(shape(mat1) .NE. shape(mat2)))  &
         stop 'NNDV_dp: shapes of input matrix 1 and input matrix 2 are not matching'
    if (any(shape(mat1) .NE. shapemask))    &
         stop 'NNDV_dp: shapes of input matrices and mask are not matching'
    !
    ! additional 2 rows and 2 cols added for checking the border cells without crashing the search agorithm
    ! so the search windows can always cover 9 cells even if it is checking a border cell
    ! buffer rows are initialized as false values within mask, i.e. they are not considered for the calculation
    ! of the criteria
    !
    ! initialize mask with default=.false.
    maske = .false.
    if (present(mask)) then
       maske(2:(size(maske,dim=1)-1_i4), 2:(size(maske,dim=2)-1_i4)) = mask
    else
       maske(2:(size(maske,dim=1)-1_i4), 2:(size(maske,dim=2)-1_i4)) = .true.
    end if
    !
    ! initialize bufferedMat1 & bufferedMat2
    bufferedMat1 = 0.0_dp
    bufferedMat2 = 0.0_dp
    bufferedMat1(2:(size(maske,dim=1)-1), 2:(size(maske,dim=2)-1)) = mat1
    bufferedMat2(2:(size(maske,dim=1)-1), 2:(size(maske,dim=2)-1)) = mat2
    !
    ! initialize NNDV
    NNDVMatrix = 0.0_dp
    !
    NNDV_dp = 0.0_dp
    do iCo = 2_i4, size(bufferedMat1, dim=2) - 1
       do iRo = 2_i4, size(bufferedMat1, dim=1) - 1
          if (.NOT. maske(iRo,iCo)) cycle
          NNDVMatrix(iRo-1_i4,iCo-1_i4) =   &
               real(    &
                  abs(  &
                     count((bufferedMat1(iRo-1:iRo+1 , iCo-1:iCo+1) - bufferedMat1(iRo,iCo) > epsilon(0.0_dp)) .AND. &
                           (maske(iRo-1:iRo+1 , iCo-1:iCo+1))) -                                                     &
                     count((bufferedMat2(iRo-1:iRo+1 , iCo-1:iCo+1) - bufferedMat2(iRo,iCo) > epsilon(0.0_dp)) .AND. &
                           (maske(iRo-1:iRo+1 , iCo-1:iCo+1)))                                                       &
                  ),    &
               dp)
          ! count - 1 to exclude referendce cell (iRo, iCo)
          validcount(iRo-1_i4,iCo-1_i4) = count(maske(iRo-1_i4:iRo+1_i4 , iCo-1_i4:iCo+1_i4)) - 1_i4
        end do
    end do
    !
    ! normalize every pixel to number of valid neighbouring cells (defined by maske) [0..8] --> [0..1]
    validmaske = (maske(2:size(maske, dim=1) - 1_i4, 2:size(maske, dim=2) - 1_i4) .and. (validcount > 0_i4))
    noValidPixels = count(validmaske)
    if (noValidPixels .GT. 0_i4) then
       NNDVMatrix = merge(NNDVMatrix/real(validcount, dp),NNDVMatrix, validmaske)
       ! average over all pixels
       NNDV_dp = 1.0_dp - sum(NNDVMatrix, mask=validmaske) / noValidPixels
       if (present(valid)) valid = .TRUE.
    ! case if maske is full of .false. or no valid neighbouring cells available for every pixel (validcount(:,:)=0)
    else
       NNDV_dp = 0.0_dp
       if (present(valid)) valid = .FALSE.
    end if
    !
  END FUNCTION NNDV_dp

  ! ----------------------------------------------------------------------------------------------------------------

  FUNCTION PD_sp(mat1, mat2, mask, valid)

    IMPLICIT NONE

    REAL(sp),    DIMENSION(:,:),                                            INTENT(IN)  :: mat1, mat2
    LOGICAL,     DIMENSION(:,:)                                 , OPTIONAL, INTENT(IN)  :: mask
    LOGICAL,                                                      OPTIONAL, INTENT(OUT) :: valid
    REAL(sp)                                                                            :: PD_sp

    INTEGER(i4)                                                                         :: iCo, iRo
    INTEGER(i4)                                                                         :: noValidPixels
    INTEGER(i4), DIMENSION(size(shape(mat1)) )                                          :: shapemask
    INTEGER(i4), DIMENSION(size(mat1, dim=1), size(mat1, dim=2))                        :: validcount
    REAL(sp),    DIMENSION(size(mat1, dim=1)+2_i4, size(mat1, dim=2)+2_i4)              :: bufferedMat1, bufferedMat2
    REAL(sp),    DIMENSION(size(mat1, dim=1), size(mat1, dim=2))                        :: PDMatrix
    LOGICAL,     DIMENSION(size(mat1, dim=1)+2_i4, size(mat1, dim=2)+2_i4)              :: maske
    LOGICAL,     DIMENSION(size(mat1, dim=1), size(mat1, dim=2))                        :: validmaske
    ! check if input has all the same dimensions
    if (present(mask)) then
       shapemask = shape(mask)
    else
       shapemask = shape(mat1)
    end if
    !
    if (any(shape(mat1) .NE. shape(mat2)))  &
         stop 'PD_sp: shapes of input matrix 1 and input matrix 2 are not matching'
    if (any(shape(mat1) .NE. shapemask))    &
         stop 'PD_sp: shapes of input matrices and mask are not matching'
    !
    ! additional 2 rows and 2 cols added for checking the border cells without crashing the search agorithm
    ! so the search windows can always cover 9 cells even if it is checking a border cell
    ! buffer rows are initialized as false values within mask, i.e. they are not considered for the calculation
    ! of the criteria
    !
    ! initialize mask with default=.false.
    maske = .false.
    if (present(mask)) then
       maske(2:(size(maske,dim=1)-1_i4), 2:(size(maske,dim=2)-1_i4)) = mask
    else
       maske(2:(size(maske,dim=1)-1_i4), 2:(size(maske,dim=2)-1_i4)) = .true.
    end if
    !
    ! initialize bufferedMat1 & bufferedMat2
    bufferedMat1 = 0.0_sp
    bufferedMat2 = 0.0_sp
    bufferedMat1(2:(size(maske,dim=1)-1), 2:(size(maske,dim=2)-1)) = mat1
    bufferedMat2(2:(size(maske,dim=1)-1), 2:(size(maske,dim=2)-1)) = mat2
    !
    ! initialize PD
    PDMatrix = 0.0_sp
    do iCo = 2_i4, size(bufferedMat1, dim=2)-1_i4
       do iRo = 2_i4, size(bufferedMat1, dim=1)-1_i4
          ! no calculation at unmasked values
          if (.NOT. maske(iRo,iCo)) cycle
          ! .NEQV. is the fortran standard for .XOR.
          ! result is written to -1 column and row because of buffer
          PDMatrix(iRo-1_i4,iCo-1_i4) = &
            real( &
              count( &
                     ( &
                        ! determine higher neighbouring values in mat1
                        (bufferedMat1(iRo-1_i4:iRo+1_i4, iCo-1_i4:iCo+1_i4) -                         &
                         bufferedMat1(iRo,iCo)                              > epsilon(0.0_sp)) .NEQV. &
                        ! determine higher neighbouring values in mat2                                &
                        (bufferedMat2(iRo-1_i4:iRo+1_i4, iCo-1_i4:iCo+1_i4) -                         &
                         bufferedMat2(iRo,iCo)                              > epsilon(0.0_sp))        &
                     ) &
                  ! exclude unmasked values
                  .and. (maske(iRo-1_i4:iRo+1_i4 , iCo-1_i4:iCo+1_i4)) &
                  ), &
            sp )
          ! count - 1 to exclude reference cell / center pixel (iRo, iCo)
          validcount(iRo-1_i4,iCo-1_i4) = count(maske(iRo-1_i4:iRo+1_i4 , iCo-1_i4:iCo+1_i4)) - 1_i4
          !
       end do
    end do
    !
    ! normalize every pixel to number of valid neighbouring cells (defined by maske) [0..8] --> [0..1]
    validmaske = (maske(2:size(maske, dim=1) - 1_i4, 2:size(maske, dim=2) - 1_i4) .and. (validcount > 0_i4))
    noValidPixels = count(validmaske)
    if (noValidPixels .GT. 0_i4) then
       PDMatrix = merge(PDMatrix/real(validcount, sp),PDMatrix, validmaske)
       ! average over all pixels
       PD_sp = 1.0_sp - sum(PDMatrix, mask=validmaske) / noValidPixels
       if (present(valid)) valid = .TRUE.
    ! case if maske is full of .false. or no valid neighbouring cells available for every pixel (validcount(:,:)=0)
    else
       PD_sp = 0.0_sp
       if (present(valid)) valid = .FALSE.
    end if
    !
  END FUNCTION PD_sp

  FUNCTION PD_dp(mat1, mat2, mask, valid)

    IMPLICIT NONE

    REAL(dp),    DIMENSION(:,:),                                            INTENT(IN)  :: mat1, mat2
    LOGICAL,     DIMENSION(:,:)                                 , OPTIONAL, INTENT(IN)  :: mask
    LOGICAL,                                                      OPTIONAL, INTENT(OUT) :: valid
    REAL(dp)                                                                            :: PD_dp

    INTEGER(i4)                                                                         :: iCo, iRo
    INTEGER(i4)                                                                         :: noValidPixels
    INTEGER(i4), DIMENSION(size(shape(mat1)) )                                          :: shapemask
    INTEGER(i4), DIMENSION(size(mat1, dim=1), size(mat1, dim=2))                        :: validcount
    REAL(dp),    DIMENSION(size(mat1, dim=1)+2_i4, size(mat1, dim=2)+2_i4)              :: bufferedMat1, bufferedMat2
    REAL(dp),    DIMENSION(size(mat1, dim=1), size(mat1, dim=2))                        :: PDMatrix
    LOGICAL,     DIMENSION(size(mat1, dim=1)+2_i4, size(mat1, dim=2)+2_i4)              :: maske
    LOGICAL,     DIMENSION(size(mat1, dim=1), size(mat1, dim=2))                        :: validmaske
    ! check if input has all the same dimensions
    if (present(mask)) then
       shapemask = shape(mask)
    else
       shapemask = shape(mat1)
    end if
    !
    if (any(shape(mat1) .NE. shape(mat2)))  &
         stop 'PD_dp: shapes of input matrix 1 and input matrix 2 are not matching'
    if (any(shape(mat1) .NE. shapemask))    &
         stop 'PD_dp: shapes of input matrices and mask are not matching'
    !
    ! additional 2 rows and 2 cols added for checking the border cells without crashing the search agorithm
    ! so the search windows can always cover 9 cells even if it is checking a border cell
    ! buffer rows are initialized as false values within mask, i.e. they are not considered for the calculation
    ! of the criteria
    !
    ! initialize mask with default=.false.
    maske = .false.
    if (present(mask)) then
       maske(2:(size(maske,dim=1)-1_i4), 2:(size(maske,dim=2)-1_i4)) = mask
    else
       maske(2:(size(maske,dim=1)-1_i4), 2:(size(maske,dim=2)-1_i4)) = .true.
    end if
    !
    ! initialize bufferedMat1 & bufferedMat2
    bufferedMat1 = 0.0_dp
    bufferedMat2 = 0.0_dp
    bufferedMat1(2:(size(maske,dim=1)-1), 2:(size(maske,dim=2)-1)) = mat1
    bufferedMat2(2:(size(maske,dim=1)-1), 2:(size(maske,dim=2)-1)) = mat2
    !
    ! initialize PD
    PDMatrix = 0.0_dp
    do iCo = 2_i4, size(bufferedMat1, dim=2)-1_i4
       do iRo = 2_i4, size(bufferedMat1, dim=1)-1_i4
          ! no calculation at unmasked values
          if (.NOT. maske(iRo,iCo)) cycle
          ! .NEQV. is the fortran standard for .XOR.
          ! result is written to -1 column and row because of buffer
          PDMatrix(iRo-1_i4,iCo-1_i4) = &
            real( &
              count( &
                     ( &
                        ! determine higher neighbouring values in mat1
                        (bufferedMat1(iRo-1_i4:iRo+1_i4, iCo-1_i4:iCo+1_i4) -                         &
                         bufferedMat1(iRo,iCo)                              > epsilon(0.0_dp)) .NEQV. &
                        ! determine higher neighbouring values in mat2                                &
                        (bufferedMat2(iRo-1_i4:iRo+1_i4, iCo-1_i4:iCo+1_i4) -                         &
                         bufferedMat2(iRo,iCo)                              > epsilon(0.0_dp))        &
                     ) &
                  ! exclude unmasked values
                  .and. (maske(iRo-1_i4:iRo+1_i4 , iCo-1_i4:iCo+1_i4)) &
                  ), &
            dp )
          ! count - 1 to exclude referendce cell (iRo, iCo)
          validcount(iRo-1_i4,iCo-1_i4) = count(maske(iRo-1_i4:iRo+1_i4 , iCo-1_i4:iCo+1_i4)) - 1_i4
          !
       end do
    end do
    !
    ! normalize every pixel to number of valid neighbouring cells (defined by maske) [0..8] --> [0..1]
    validmaske = (maske(2:size(maske, dim=1) - 1_i4, 2:size(maske, dim=2) - 1_i4) .and. (validcount > 0_i4))
    noValidPixels = count(validmaske)
    if (noValidPixels .GT. 0_i4) then
       PDMatrix = merge(PDMatrix/real(validcount, dp),PDMatrix, validmaske)
       ! average over all pixels
       PD_dp = 1.0_dp - sum(PDMatrix, mask=validmaske) / noValidPixels
       if (present(valid)) valid = .TRUE.
    ! case if maske is full of .false. or no valid neighbouring cells available for every pixel (validcount(:,:)=0)
    else
       PD_dp = 0.0_dp
       if (present(valid)) valid = .FALSE.
    end if
    !
  END FUNCTION PD_dp
  !
END MODULE mo_spatialsimilarity


module bar
    implicit none
    public
    integer, allocatable :: global_var
    !$acc declare present(global_var)
contains
end module

subroutine foo(sub_data, n)
    !$acc routine vector
    use bar
    implicit none

    integer i, n
    integer sub_data(n)
    print*, 'GPU thread reads a variable from module \"bar\"'
    !$acc loop
    do i = 1, n
        sub_data(i) = sub_data(i) + 1 + global_var
    end do
end subroutine

program main
    !$acc routine(foo) vector
    use bar
    implicit none

    integer y, w, h
    integer, allocatable :: all_data(:, :)

    w = 6
    h = 3
    allocate(all_data(w, h))
    all_data = 0

    global_var = 123

    !$acc parallel loop gang
    do y = 1, h
        all_data(:, y) = y
        call foo(all_data(:, y), w)
    end do

    print*, all_data
    deallocate(all_data)
end program main


text:
! { dg-do compile }
! Fixes PR37787 where the arithmetic overflow was not detected and an ICE ensued.
!
! Contributed by Tobias Burnus <burnus@gcc.gnu.org>
!
program bug
  implicit none
   integer(1) :: a(2) = (/ Z'FF', Z'FF' /) ! { dg-error "Arithmetic overflow" }
   print*, a
end program bug

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
       DOUBLE PRECISION FUNCTION getMass(Ptype)
c###############################################
c                                             ##
c         Masses of Particles in GeV          ##
c                                             ##
c###############################################
         IMPLICIT NONE
         INTEGER          Ptype

         SELECT CASE(Ptype)
            CASE( 0)                                       !________ B O S O N S
                   getMass = 0d0             ! Gamma
            CASE(-1)
                   getMass = 91.1867d0       ! Z
            CASE(-2)
                   getMass = 120d0           ! H
            CASE(-3)
                   getMass = 80.4514958d0    ! W^{+}
            CASE(-4)
                   getMass = 80.4514958d0    ! W^{-}
                                                           !________ L E P T O N T S
            CASE( 1)
                   getMass = 0d0             ! neutrino_electron
            CASE( 2)
                   getMass = 0.51099907d-3   ! electron
            CASE( 3)
                   getMass = 0d0             ! neutrino_muon
            CASE( 4)
                   getMass = 105.658389d-3   ! miuon
            CASE( 5)
                   getMass = 0d0             ! neutrino_tau
            CASE( 6)
                   getMass = 1.77705d0       ! tau

            CASE( 7)                                       !________ Q U A R K S
                   getMass = 62.0d-3         ! up
            CASE( 8)
                   getMass = 83.0d-3         ! down
            CASE( 9)
                   getMass = 1.5d0           ! c
            CASE(10)
                   getMass = 215.0d-3        ! s
            CASE(11)
                   getMass = 173.8d0         ! t
            CASE(12)
                   getMass = 4.7d0           ! b
            CASE default
               CALL MessageRoutine(1,.true.,"getMass")
               STOP
         END SELECT

       RETURN
       END

       DOUBLE PRECISION FUNCTION getCharge(Ptype)
c############################################################
c                                                          ##
c   Charges for Particles (in the electron charge unit!!)  ##
c                                                          ##
c############################################################
         IMPLICIT NONE
         INTEGER          Ptype

         SELECT CASE(Ptype)
            CASE( 0)         ! Gamma                       !________ B O S O N S
                   getCharge =  0d0
            CASE(-1)
                   getCharge =  0d0         ! Z
            CASE(-2)
                   getCharge =  0d0         ! H
            CASE(-3)
                   getCharge = +1d0         ! W^{+}
            CASE(-4)
                   getCharge = -1d0         ! W^{-}
                                                           !________ L E P T O N T S
            CASE( 1)
                   getCharge =  0d0         ! neutrino_electron
            CASE( 2)
                   getCharge = -1d0         ! electron
            CASE( 3)
                   getCharge =  0d0         ! neutrino_muon
            CASE( 4)
                   getCharge = -1d0         ! miuon
            CASE( 5)
                   getCharge =  0d0         ! neutrino_tau
            CASE( 6)
                   getCharge = -1d0         ! tau
                                                            !________ Q U A R K S
            CASE( 7)
                   getCharge = +2d0/3d0     ! up
            CASE( 8)
                   getCharge = -1d0/3d0     ! down
            CASE( 9)
                   getCharge = +2d0/3d0     ! c
            CASE(10)
                   getCharge = -1d0/3d0     ! s
            CASE(11)
                   getCharge = +2d0/3d0     ! t
            CASE(12)
                   getCharge = -1d0/3d0     ! b
            CASE default
               CALL MessageRoutine(1,.true.,"getCharge")
               STOP
         END SELECT

       RETURN
       END

       DOUBLE PRECISION FUNCTION getWIS(Ptype)
c##################################################
c                                                ##
c   Weak IsoSpin (I3i) of Leptons and Quarks ... ##
c                                                ##
c##################################################
         IMPLICIT NONE
         INTEGER          Ptype

         SELECT CASE(Ptype)
                                                            !________ L E P T O N T S
            CASE( 1)
                   getWIS = +1d0/2d0      ! neutrino_electron
            CASE( 2)
                   getWIS = -1d0/2d0      ! electron
            CASE( 3)
                   getWIS = +1d0/2d0      ! neutrino_muon
            CASE( 4)
                   getWIS = -1d0/2d0      ! miuon
            CASE( 5)
                   getWIS = +1d0/2d0      ! neutrino_tau
            CASE( 6)
                   getWIS = -1d0/2d0       ! tau
                                                            !________ Q U A R K S
            CASE( 7)
                   getWIS = +1d0/2d0      ! up
            CASE( 8)
                   getWIS = -1d0/2d0      ! down
            CASE( 9)
                   getWIS = +1d0/2d0      ! c
            CASE(10)
                   getWIS = -1d0/2d0      ! s
            CASE(11)
                   getWIS = +1d0/2d0      ! t
            CASE(12)
                   getWIS = -1d0/2d0      ! b
            CASE default
               CALL MessageRoutine(1,.true.,"getWIS")
               STOP
         END SELECT

       RETURN
       END

       SUBROUTINE ConvertPNS(upPart,dnPart,s2nIn)
c######################################################
c    Converts Particle Numbering Scheme from S2N to   #
c    Monte Carlo                                      #
c                                                     #
c    Input : upPart, dnPart                           #
c    Output: s2nIn                                    #
c    s2nIn is a channel identification in s2n !!!     #
c######################################################
         IMPLICIT NONE
         INTEGER s2nIn,sxv1,sxv2,upPart,dnPart,dcPart
         INCLUDE "MC_Declare.h"

         sxv1 = dnPart-upPart
         sxv2 = dnPart-upPart+1-2*mod(upPart,2)          ! check of doublets...
         IF (sxv1.NE.0 .AND. (decPart.EQ.-1 .OR. decPart.EQ.-2)) THEN
           CALL MessageRoutine(1,.true.,"ConvertPNS")
           STOP
         ENDIF
         IF (sxv2.NE.0 .AND. (decPart.EQ.-3 .OR. decPart.EQ.-4)) THEN
           CALL MessageRoutine(1,.true.,"ConvertPNS")
           STOP
         ENDIF
         IF ((upPart.GE.7.AND.upPart.LE.12) .AND. (dnPart.GE.7.AND.dnPart.LE.12)) THEN
            s2nIn = Min(upPart,dnPart)-3
         ELSEIF ((upPart.GE.1.AND.upPart.LE.6) .AND. (dnPart.GE.1.AND.dnPart.LE.6)) THEN
            s2nIn = Max(upPart,dnPart)/2
         ELSE
            CALL MessageRoutine(1,.true.,"CoonvertPNS")
            STOP
         ENDIF

       RETURN
       END

      SUBROUTINE  MessageRoutine(type,error,procedureName)

        IMPLICIT NONE
        CHARACTER procedureName*(*),genType*30,softEvent*3,hardEvent*3
        LOGICAL   error
        INTEGER   type
        INCLUDE   "MC_Declare.h"

        IF (type.EQ.1) THEN       ! notation of particles ....

           WRITE(mcLUN,*)""
           WRITE(mcLUN,*)"     ERROR IN  ",procedureName," -- USE CORRECT VALUE FOR `Ptype':"
           WRITE(mcLUN,*)""

           IF (error.AND.type.EQ.1) THEN
             WRITE(mcLUN,*)""
             WRITE(mcLUN,'(3X,A65)') "#################__PARTICLE__NUMBERING__SCHEME__#################"
             WRITE(mcLUN,'(3X,A65)') "#                                                               #"
             WRITE(mcLUN,'(3X,A65)') "# _B_O_S_O_N_S_:                                                #"
             WRITE(mcLUN,'(3X,A65)') "#   gamma = 0  Z = -1  H = -2  W^{+} = -3  W^{-} = -4           #"
             WRITE(mcLUN,'(3X,A65)') "#                                                               #"
             WRITE(mcLUN,'(3X,A65)') "# _L_E_P_T_O_N_S_:                                              #"
             WRITE(mcLUN,'(3X,A65)') "#   neutrino_electron = 1  neutrino_muon = 3  neutrino_tau = 5  #"
             WRITE(mcLUN,'(3X,A65)') "#   electron          = 2  muon          = 4  tau          = 6  #"
             WRITE(mcLUN,'(3X,A65)') "#                                                               #"
             WRITE(mcLUN,'(3X,A65)') "# _Q_U_A_R_K_S_:                                                #"
             WRITE(mcLUN,'(3X,A65)') "#   up   = 7  c = 9   t = 11                                    #"
             WRITE(mcLUN,'(3X,A65)') "#   down = 8  s = 10  b = 12                                    #"
             WRITE(mcLUN,'(3X,A65)') "#################################################################"
             WRITE(mcLUN,'(3X,A65)')""
           ENDIF

        ELSEIF (type.EQ.2.AND.error) THEN           ! random number generator..
           IF (ranGenType .EQ. 1) THEN
              genType = "__RANLUX__"
           ELSEIF (ranGenType .EQ. 2) THEN
              genType = "__RANMAR__"
           ELSEIF (ranGenType .EQ. 3) THEN
              genType = "__MERSENNE TWISTER__"
           ELSE
              WRITE(mcLUN,*)""
              WRITE(mcLUN,*)"  ERROR IN  ",procedureName,"  :"
              WRITE(mcLUN,*)"        NO RANDOM GNERATOR IS INITIALIZED:"
              WRITE(mcLUN,*)"            RanGenType = 1 for 'ranlux'"
              WRITE(mcLUN,*)"            RanGenType = 2 for 'ranmar'"
              WRITE(mcLUN,*)"            RanGenType = 3 for 'mersenne twister'"
              WRITE(mcLUN,*)""
              STOP
           ENDIF
           WRITE(mcLUN,*) ""
           WRITE(mcLUN,'(3X,A42)') "##########################################"
           WRITE(mcLUN,'(3X,A42)') "#                                        #"
           WRITE(mcLUN,'(3X,A42)') "# USED RANDOM NUMBER GENERATOR :         #"
           WRITE(mcLUN,'(3X,A42)') "#                                        #"
           WRITE(mcLUN,'(3X,A10,A30,A2)') "#         ",genType,"#"
           WRITE(mcLUN,'(3X,A42)') "##########################################"
           WRITE(mcLUN,*) ""

        ELSEIF (type.EQ.3.AND.error) THEN           ! Enevt type messeger ...
           IF(.NOT.softevnt.AND..NOT.hardEvnt) THEN
             WRITE(mcLUN,*) ""
             WRITE(mcLUN,*) "ERROR IN ",procedureName,"   :"
             WRITE(mcLUN,*) "      NO EVENTS WILL BE GENERATED :"
             WRITE(mcLUN,*) "             Hard Event =",hardEvnt
             WRITE(mcLUN,*) "             Soft Event =",softEvnt
             WRITE(mcLUN,*) ""
             STOP
           ENDIF
           IF (hardEvnt)      hardEvent="ON"
           IF (.NOT.hardEvnt) hardEvent="OFF"
           IF (softEvnt)      softEvent="ON"
           IF (.NOT.softEvnt) softEvent="OFF"
           WRITE(mcLUN,*) ""
           WRITE(mcLUN,*) "MESSAGE FROM ",procedureName,"   :"
           WRITE(mcLUN,*) "        TYPE OF GENERATED EVENTS :"
           WRITE(mcLUN,*) "               Hard Event = ",hardEvent
           WRITE(mcLUN,*) "               Soft Event = ",softEvent
           WRITE(mcLUN,*) ""

        ELSEIF (type.EQ.4.AND.error) THEN           ! Check event weight sign on negativeness !!!
           WRITE(mcLUN,*) " "
           WRITE(mcLUN,*) " ERROR IN ",procedurename,"  :"
           WRITE(mcLUN,*) "       PRODUCED NEGATIVE WEIGHT !!!...k0Min IS TOO SMALL !!!"
           WRITE(mcLUN,*) ""
           STOP

        ELSEIF (type.EQ.5.AND.error) THEN           ! ElectroWeak Correction ON|OFF !!!
           IF (EWC .EQ. 1d0) THEN
              WRITE(mcLUN,*) ""
              WRITE(mcLUN,'(3X,A40)') "########################################"
              WRITE(mcLUN,'(3X,A40)') "#                                      #"
              WRITE(mcLUN,'(3X,A40)') "# MESSAGE FROM MC_Init :               #"
              WRITE(mcLUN,'(3X,A40)') "#         EW CORRECTIONS -- SWITCH ON  #"
              WRITE(mcLUN,'(3X,A40)') "#                                      #"
              WRITE(mcLUN,'(3X,A40)') "########################################"
              WRITE(mcLUN,*) ""
           ELSEIF (EWC.EQ.0d0) THEN
              WRITE(mcLUN,*) ""
              WRITE(mcLUN,'(3X,A40)') "########################################"
              WRITE(mcLUN,'(3X,A40)') "#                                      #"
              WRITE(mcLUN,'(3X,A40)') "# MESSAGE FROM MC_Init :               #"
              WRITE(mcLUN,'(3X,A40)') "#         EW CORRECTIONS -- SWITCH OFF #"
              WRITE(mcLUN,'(3X,A40)') "#                                      #"
              WRITE(mcLUN,'(3X,A40)') "########################################"
              WRITE(mcLUN,*) ""
           ELSE
              WRITE(mcLUN,*) ""
              WRITE(mcLUN,'(3X,A38)') "######################################"
              WRITE(mcLUN,'(3X,A38)') "#                                    #"
              WRITE(mcLUN,'(3X,A38)') "#   ERROR IN MC_Init :               #"
              WRITE(mcLUN,'(3X,A38)') "#         EWC FLAG MUST BE :         #"
              WRITE(mcLUN,'(3X,A38)') "#          1 -- EW correction ON     #"
              WRITE(mcLUN,'(3X,A38)') "#          0 -- EW correction OFF    #"
              WRITE(mcLUN,'(3X,A38)') "#                                    #"
              WRITE(mcLUN,'(3X,A38)') "######################################"
              WRITE(mcLUN,*) ""
              STOP
           ENDIF
        ELSEIF (type.EQ.6.AND.error) THEN           ! Photon Energy Range Checks !!!

           IF ( k0Max.GT.k0MaxKin ) THEN
              PRINT*, ""
              PRINT*, "ERROR IN MakePhaseSpace_3body : "
              PRINT*, "      INCORRECT VALUE FOR k0Max "
              PRINT*, "             k0Max       = ",k0Max
              PRINT*, "             k0Max MUST <= ",k0MaxKin
              PRINT*, ""
              STOP
           ELSEIF ( omega.GT.k0MaxKin ) THEN
              PRINT*, ""
              PRINT*, "ERROR IN MakePhaseSpace_3body : "
              PRINT*, "      INCORRECT VALUE FOR omega "
              PRINT*, "             omega    = ",omega
              PRINT*, "          omega MUST <= ",k0MaxKin
              PRINT*, ""
              STOP
          ELSEIF ( k0Max.LT.omega ) THEN
              PRINT*, ""
              PRINT*, "ERROR IN MakePhaseSpace_3body : "
              PRINT*, "      INCORRECT VALUE FOR omega  "
              PRINT*, "           omega    = ",omega
              PRINT*, "        omega MUST <= ",k0Max
              PRINT*, ""
              STOP
          ENDIF

        ENDIF
      RETURN
      END


text:
!
! Test_TauCoeff
!
! Program to test the TauCoeff definition and I/O routines.
!
!
! FILES ACCESSED:
!       - netCDF format TauCoeff data file, both input and output
!       - Binary format TauCoeff data file, both input and output
!
! SIDE EFFECTS:
!       The test output files are overwritten if they already exist.
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, CIMSS/SSEC 02-Jan-2003
!                       paul.vandelst@ssec.wisc.edu
!

PROGRAM Test_TauCoeff

  ! -----------------
  ! Environment setup
  ! -----------------
  ! Module usage
  USE Type_Kinds
  USE Message_Handler
  USE TauCoeff_Define
  USE TauCoeff_Binary_IO
  USE TauCoeff_netCDF_IO
  ! Disable implicit typing
  IMPLICIT NONE

  ! ----------
  ! Parameters
  ! ----------
  CHARACTER(*), PARAMETER :: PROGRAM_NAME = 'Test_TauCoeff'
  CHARACTER(*), PARAMETER :: PROGRAM_RCS_ID = &
  CHARACTER(*), PARAMETER ::  NC_FILENAME = 'Test.TauCoeff.nc'
  CHARACTER(*), PARAMETER :: BIN_FILENAME = 'Test.TauCoeff.bin'
  INTEGER, PARAMETER :: SET = 1
  INTEGER, PARAMETER :: MAX_N_LOOPS  = 5000
  INTEGER, PARAMETER :: INFO_N_LOOPS = 500

  ! ---------
  ! Variables
  ! ---------
  CHARACTER(256) :: Message
  INTEGER :: Error_Status
  INTEGER :: n_Orders
  INTEGER :: n_Predictors, i
  INTEGER :: n_Absorbers,  j
  INTEGER :: n_Channels,   l
  INTEGER :: n
  TYPE(TauCoeff_type) :: TauCoeff1
  TYPE(TauCoeff_type) :: TauCoeff2
  TYPE(TauCoeff_type) :: TauCoeff3

  ! Program header
  ! --------------
  CALL Program_Message( PROGRAM_NAME, &
                        'Program to test TauCoeff definition and I/O routines', &
                        '$Revision$' )

  ! Fill a structure with dummy values
  ! ----------------------------------
  ! Define some typical dimensions
  n_Orders     = 10
  n_Predictors = 5
  n_Absorbers  = 3
  n_Channels   = 100

  ! Allocate a TauCoeff structure
  Error_Status = Allocate_TauCoeff( n_Orders, &
                                    n_Predictors, &
                                    n_Absorbers, &
                                    n_Channels, &
                                    TauCoeff1 )
  IF ( Error_Status /= SUCCESS ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'Error allocating TauCoeff1 structre.', &
                          Error_Status )
    STOP
  END IF

  ! Fill up the structure with stuff
  TauCoeff1%n_Sensors = 1
  TauCoeff1%Sensor_Descriptor = 'sensor_platform'
  TauCoeff1%NCEP_Sensor_ID   = 1
  TauCoeff1%WMO_Satellite_ID = 2
  TauCoeff1%WMO_Sensor_ID    = 3
  TauCoeff1%Sensor_Channel   = (/(l,l=1,n_Channels)/)
  TauCoeff1%Absorber_ID = (/(j,j=1,n_Absorbers)/)
  TauCoeff1%Alpha    = 4.0_Double
  TauCoeff1%Alpha_C1 = 5.0_Double
  TauCoeff1%Alpha_C2 = 6.0_Double
  TauCoeff1%Order_Index = 5
  DO l = 1, n_Channels
    DO j = 1, n_Absorbers
      TauCoeff1%Predictor_Index(:,j,l) = (/(i,i=0,n_Predictors)/)
    END DO
  END DO
  TauCoeff1%C = -99.0_Double

  ! Write and read a netCDF TauCoeff file
  ! -------------------------------------
  WRITE( *, '( /5x, "Testing TauCoeff netCDF I/O functions ..." )' )

  ! Write the netCDF data file
  Error_Status = Write_TauCoeff_netCDF( NC_FILENAME, &
                                        TauCoeff1, &
                                        Title        ='This is the title attribute', &
                                        History      ='This is the history attribute', &
                                        Sensor_Name  ='This is the sensor_name attribute',   &
                                        Platform_Name='This is the platform_name attribute', &
                                        Comment      ='This is the comment attribute',&
                                        ID_Tag       ='This is the id_tag attribute' )
  IF ( Error_Status /= SUCCESS ) THEN
    Error_Status = WARNING
    Message = 'Error writing the netCDF TauCoeff file '//TRIM(NC_FILENAME)
  ELSE
    Error_Status = INFORMATION
    Message = TRIM(NC_FILENAME)//' netCDF file write was successful'
  END IF
  CALL Display_Message( PROGRAM_NAME, &
                        TRIM(Message), &
                        Error_Status )

  ! Read the netCDF data file
  Error_Status = Read_TauCoeff_netCDF( NC_FILENAME, &
                                       TauCoeff2 )
  IF ( Error_Status /= SUCCESS ) THEN
    Error_Status = WARNING
    Message = 'Error reading the netCDF TauCoeff file '//TRIM(NC_FILENAME)
  ELSE
    Error_Status = INFORMATION
    Message = TRIM(NC_FILENAME)//' netCDF file read was successful'
  END IF
  CALL Display_Message( PROGRAM_NAME, &
                        TRIM(Message), &
                        Error_Status )

  ! Compare the structures
  Error_Status = Equal_TauCoeff( TauCoeff2, TauCoeff1 )
  IF ( Error_Status /= SUCCESS ) THEN
    Error_Status = WARNING
    Message = 'NetCDF I/O structures are different'
  ELSE
    Error_Status = INFORMATION
    Message = 'NetCDF I/O structures are equal'
  END IF

  CALL Display_Message( PROGRAM_NAME, &
                        TRIM(Message), &
                        Error_Status )

  WRITE( *,'(//5x,"Press <ENTER> to test Binary I/O functions...")' )
  READ( *,* )

  ! Write and read a Binary TauCoeff file
  ! -------------------------------------
  WRITE( *, '( /5x, "Testing TauCoeff Binary Write/Read functions ..." )' )

  ! Write the Binary data file
  Error_Status = Write_TauCoeff_Binary( BIN_FILENAME, &
                                        TauCoeff1 )
  IF ( Error_Status /= SUCCESS ) THEN
    Error_Status = WARNING
    Message = 'Error writing the Binary TauCoeff file '//TRIM(BIN_FILENAME)
  ELSE
    Error_Status = INFORMATION
    Message = TRIM(BIN_FILENAME)//' Binary file write was successful'
  END IF
  CALL Display_Message( PROGRAM_NAME, &
                        TRIM(Message), &
                        Error_Status )

  ! Read the Binary data file
  Error_Status = Read_TauCoeff_Binary( BIN_FILENAME, &
                                       TauCoeff3 )
  IF ( Error_Status /= SUCCESS ) THEN
    Error_Status = WARNING
    Message = 'Error reading the Binary TauCoeff file '//TRIM(BIN_FILENAME)
  ELSE
    Error_Status = INFORMATION
    Message = TRIM(BIN_FILENAME)//' Binary file read was successful'
  END IF
  CALL Display_Message( PROGRAM_NAME, &
                        TRIM(Message), &
                        Error_Status )

  ! Compare the structures
  Error_Status = Equal_TauCoeff( TauCoeff3, TauCoeff1 )
  IF ( Error_Status /= SUCCESS ) THEN
    Error_Status = WARNING
    Message = 'Binary I/O structures are different'
  ELSE
    Error_Status = INFORMATION
    Message = 'Binary I/O structures are equal'
  END IF
  CALL Display_Message( PROGRAM_NAME, &
                        TRIM(Message), &
                        Error_Status )

  WRITE( *,'(//5x,"Press <ENTER> to test for netCDF reader memory leaks...")' )
  READ( *,* )

  ! Loop for memory leak test
  ! -------------------------
  ! Test the netCDF reader for memory leaks
  WRITE( *,'(/5x,"Looping for netCDF read memory leak test ...")' )
  DO n = 1, MAX_N_LOOPS
    Error_Status = Read_TauCoeff_netCDF( NC_FILENAME, &
                                         TauCoeff1, &
                                         Quiet=SET )
    IF ( MOD( n, INFO_N_LOOPS ) == 0 .AND. n /= MAX_N_LOOPS ) THEN
      WRITE( Message, '( "Completed loop #",i0," of ", i0)' ) n, MAX_N_LOOPS
      CALL Display_Message( PROGRAM_NAME, &
                            TRIM(Message), &
                            INFORMATION )
    END IF
  END DO
  WRITE( *,'(//5x,"Press <ENTER> to test for Binary reader memory leaks...")' )
  READ( *,* )

  ! Test the Binary reader for memory leaks
  WRITE( *,'(/5x,"Looping for Binary read memory leak test ...")' )
  DO n = 1, MAX_N_LOOPS
    Error_Status = Read_TauCoeff_Binary( BIN_FILENAME, &
                                         TauCoeff2, &
                                         Quiet=SET )
    IF ( MOD( n, INFO_N_LOOPS ) == 0 .AND. n /= MAX_N_LOOPS ) THEN
      WRITE( Message, '( "Completed loop #",i0," of ", i0)' ) n, MAX_N_LOOPS
      CALL Display_Message( PROGRAM_NAME, &
                            TRIM(Message), &
                            INFORMATION )
    END IF
  END DO
  WRITE( *,'(//5x,"Press <ENTER> to test for structure assign function memory leaks...")' )
  READ( *,* )

  ! Test the Assign function for memory leaks
  WRITE( *,'(/5x,"Looping for structure copy memory leak test ...")' )
  DO n = 1, MAX_N_LOOPS
    Error_Status = Assign_TauCoeff( TauCoeff3, TauCoeff2 )
    IF ( MOD( n, INFO_N_LOOPS ) == 0 .AND. n /= MAX_N_LOOPS ) THEN
      WRITE( Message, '( "Completed loop #",i0," of ", i0)' ) n, MAX_N_LOOPS
      CALL Display_Message( PROGRAM_NAME, &
                            TRIM(Message), &
                            INFORMATION )
    END IF
  END DO

  ! Destroy the structures
  ! ----------------------
  Error_Status = Destroy_TauCoeff( TauCoeff1 )
  IF ( Error_Status /= SUCCESS ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'Error destroying TauCoeff1 structure.', &
                          WARNING )
  END IF

  Error_Status = Destroy_TauCoeff( TauCoeff2 )
  IF ( Error_Status /= SUCCESS ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'Error destroying TauCoeff2 structure.', &
                          WARNING )
  END IF

  Error_Status = Destroy_TauCoeff( TauCoeff3 )
  IF ( Error_Status /= SUCCESS ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'Error destroying TauCoeff3 structure.', &
                          WARNING )
  END IF

END PROGRAM Test_TauCoeff

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
C Copyright(C) 2009 Sandia Corporation. Under the terms of Contract
C DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
C certain rights in this software.
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C     * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C     * Redistributions in binary form must reproduce the above
C       copyright notice, this list of conditions and the following
C       disclaimer in the documentation and/or other materials provided
C       with the distribution.
C     * Neither the name of Sandia Corporation nor the names of its
C       contributors may be used to endorse or promote products derived
C       from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

C $Log: chkerr.f,v $
C Revision 1.3  2009/03/25 12:36:42  gdsjaar
C Add copyright and license notice to all files.
C Permission to assert copyright has been granted; blot is now open source, BSD
C
C Revision 1.2  1997/11/11 14:55:53  gdsjaar
C Added 'external blkdat' to main program to ensure that the block data
C gets linked into the executable. Wasn't happening on dec alpha
C systems.
C
C Removed unreachable lines in several routines
C
C Fixed variable name spelling in contor.f
C
C Unsplit strings that were split across lines
C
C Removed old error variables left over from exodusIIv1
C
C Upped version number
C
C Revision 1.1  1994/04/07 19:55:29  gdsjaar
C Initial checkin of ACCESS/graphics/blotII2
C
c Revision 1.3  1993/07/28  15:58:32  gdsjaar
c Fixed junky looking output
c
c Revision 1.2  1993/07/27  20:27:18  gdsjaar
c Added error checking, set to non-verbose exodus II file opening,
c cleaned up dbiv0 routine, removed some unused variables.
c
c Revision 1.1  1993/07/27  19:35:17  gdsjaar
c Added file for error checking in exodus II calls.
c
c Revision 1.3  1993/06/24  15:27:20  gdsjaar
c Added more error flags and messages.
c
c Revision 1.2  1993/03/03  17:44:44  gdsjaar
c Fixed problem with assumed length strings.
c
c Revision 1.1  1992/06/08  22:23:07  gdsjaar
c New routine to check error status of exo2 library calls
c
C=======================================================================
      subroutine chkerr (routine, caller, ierr)
C=======================================================================
      character*(*) routine, caller
      include 'exodusII.inc'
      character*80 string, path
      character*40 errors(19)
      integer ierr
      integer idebug
      data errors /
     *  'Not a netcdf id',
     *  'Too many netcdfs open',
     *  'netcdf file exists && NCNOCLOB',
     *  'Invalid Argument',
     *  'Write to read only',
     *  'Operation not allowed in data mode',
     *  'Operation not allowed in define mode',
     *  'Coordinates out of Domain',
     *  'MAXNCDIMS exceeded',
     *  'String match to name in use',
     *  'Attribute not found',
     *  'MAXNCATTRS exceeded',
     *  'Not a netcdf data type',
     *  'Invalid dimension id',
     *  'NCUNLIMITED in the wrong index',
     *  'MAXNCVARS exceeded',
     *  'Variable not found',
     *  'Action prohibited on NCGLOBAL varid',
     *  'Not a netcdf file'/

      idebug = 1
      write (path, '(A,A2,A)') CALLER, '->', ROUTINE
      lp = lenstr(path)
      if (idebug .gt. 0 .and. ierr .gt. 0)
     *  call prterr ('CMDSPEC', PATH(:LP))

      if (ierr .eq. 0) then
         continue
      else if (ierr .eq. EXBFID) then
         write (string, '(A,A)')
     $        'The EXODUS file ID is incorrect in ', path(:lp)
         call prterr ('FATAL', string)
         stop 'EXBFID'
      else if (ierr .eq. EXWARN) then
         write (string, '(A,A)')
     $        'A non-fatal error occured in ', path(:lp)
         call prterr ('WARNING', string)
      else
         if (ierr .le. 19) then
           l = lenstr(errors(ierr))
           write (string, '(A,A,A)') errors(ierr)(:l), ' ', path(:lp)
         else
           write (string, '(A,A)')
     *       'Unknown Error ', path(:lp)
         end if
         call prterr ('WARNING', string)
      end if

      return
      end

max_stars_repo_name: milljm/seacas

max_issues_repo_name: milljm/seacas

max_forks_repo_name: milljm/seacas