text:
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

MODULE m_wann_orbcomp
  !*************************************************************************
  !     Compute matrix elements of the first 16 angular momentum projectors
  !     within the MT-spheres. This is done separately for each atom type.
  !
  !     These are:
  !     1: |s><s|
  !     2: |p_x><p_x|
  !     3: |p_y><p_y|
  !     4: |p_z><p_z|
  !     5: |d_xy><d_xy|
  !     6: |d_yz><d_yz|
  !     7: |d_zx><d_zx|
  !     8: |d_x2-y2><d_x2-y2|
  !     9: |d_z2><dz2|
  !     10: |f><f|
  !     11: |f><f|
  !     12: |f><f|
  !     13: |f><f|
  !     14: |f><f|
  !     15: |f><f|
  !     16: |f><f|
  !
  !     Based on code in orb_comp by Yury Koroteev
  !
  !     Frank Freimuth
  !***********************************************************************
CONTAINS
  SUBROUTINE wann_orbcomp(atoms,usdus,jspin,&
       acof,bcof,ccof,&
       num_orbs,&
       orbs,&
       l_oc_f,&
       mmn)
    USE m_types
    IMPLICIT NONE
    !     .. scalar arguments ..
    TYPE(t_atoms),INTENT(in) :: atoms
    TYPE(t_usdus),INTENT(IN) :: usdus
    INTEGER,INTENT(IN)       :: jspin
    !     .. array arguments ..
    COMPLEX, INTENT (in)  :: ccof(-atoms%llod:,:,:,:) !ccof(-llod:llod,noccbd,nlod,natd)
    COMPLEX, INTENT (in)  :: acof(:,0:,:)!acof(noccbd,0:lmd,natd)
    COMPLEX, INTENT (in)  :: bcof(:,0:,:)!bcof(noccbd,0:lmd,natd)
    COMPLEX, INTENT (inout) :: mmn(:,:,:,:)!mmn(16,noccbd,noccbd,ntype)

    INTEGER, INTENT (in)  :: num_orbs
    INTEGER, INTENT (in)  :: orbs(:)
    LOGICAL, INTENT (in)  :: l_oc_f

    !    ..Local Scalars

    COMPLEX  ca00(SIZE(acof,1))
    COMPLEX  ca01(SIZE(acof,1))
    COMPLEX  ca02(SIZE(acof,1))
    COMPLEX  ca03(SIZE(acof,1))
    COMPLEX  ca04(SIZE(acof,1))
    COMPLEX  ca05(SIZE(acof,1))
    COMPLEX  ca06(SIZE(acof,1))
    COMPLEX  ca07(SIZE(acof,1))
    COMPLEX  ca08(SIZE(acof,1))
    COMPLEX  ca09(SIZE(acof,1))
    COMPLEX  ca10(SIZE(acof,1))
    COMPLEX  ca11(SIZE(acof,1))
    COMPLEX  ca12(SIZE(acof,1))
    COMPLEX  ca13(SIZE(acof,1))
    COMPLEX  ca14(SIZE(acof,1))
    COMPLEX  ca15(SIZE(acof,1))

    COMPLEX  cb00(SIZE(acof,1))
    COMPLEX  cb01(SIZE(acof,1))
    COMPLEX  cb02(SIZE(acof,1))
    COMPLEX  cb03(SIZE(acof,1))
    COMPLEX  cb04(SIZE(acof,1))
    COMPLEX  cb05(SIZE(acof,1))
    COMPLEX  cb06(SIZE(acof,1))
    COMPLEX  cb07(SIZE(acof,1))
    COMPLEX  cb08(SIZE(acof,1))
    COMPLEX  cb09(SIZE(acof,1))
    COMPLEX  cb10(SIZE(acof,1))
    COMPLEX  cb11(SIZE(acof,1))
    COMPLEX  cb12(SIZE(acof,1))
    COMPLEX  cb13(SIZE(acof,1))
    COMPLEX  cb14(SIZE(acof,1))
    COMPLEX  cb15(SIZE(acof,1))

    COMPLEX  cc00(SIZE(acof,1))
    COMPLEX  cc01(SIZE(acof,1))
    COMPLEX  cc02(SIZE(acof,1))
    COMPLEX  cc03(SIZE(acof,1))
    COMPLEX  cc04(SIZE(acof,1))
    COMPLEX  cc05(SIZE(acof,1))
    COMPLEX  cc06(SIZE(acof,1))
    COMPLEX  cc07(SIZE(acof,1))
    COMPLEX  cc08(SIZE(acof,1))
    COMPLEX  cc09(SIZE(acof,1))
    COMPLEX  cc10(SIZE(acof,1))
    COMPLEX  cc11(SIZE(acof,1))
    COMPLEX  cc12(SIZE(acof,1))
    COMPLEX  cc13(SIZE(acof,1))
    COMPLEX  cc14(SIZE(acof,1))
    COMPLEX  cc15(SIZE(acof,1))

    INTEGER ::  lo,orbi,rep
    REAL    :: ddn0,ddn1,ddn2,ddn3
    INTEGER :: mt,ityp,n,m,mt1,l
    REAL    :: h,g,c3,c5
    LOGICAL :: l_doit
    !      ..
    !      ..
    !      Intrinsic Function
    INTRINSIC sqrt,conjg
    !
    DATA h/0.50/ g/0.0625/
    !****************************************************
    c5 = SQRT(5.0)
    c3 = SQRT(3.0)
    !
    mt=0
    DO  ityp = 1,atoms%ntype
       DO rep=1,atoms%neq(ityp)
          mt=mt+1

          l_doit=.FALSE.
          DO l=1,num_orbs
             IF(orbs(l)==mt)THEN
                l_doit=.TRUE.
                orbi=l
             ENDIF
          ENDDO
          IF(.NOT.l_doit)CYCLE

          ddn0 = usdus%ddn(0,ityp,jspin)
          ddn1 = usdus%ddn(1,ityp,jspin)
          ddn2 = usdus%ddn(2,ityp,jspin)
          ddn3 = usdus%ddn(3,ityp,jspin)
          DO n=1,SIZE(acof,1)
             !
             !acof
             !         s-states
             ca00(n)=acof(n,0,mt)
             !        p-states
             ca01(n) = acof(n,1,mt) - acof(n,3,mt)
             ca02(n) = acof(n,1,mt) + acof(n,3,mt)
             ca03(n) = acof(n,2,mt)
             !        d-states
             ca04(n) = acof(n,4,mt) - acof(n,8,mt)
             ca05(n) = acof(n,5,mt) + acof(n,7,mt)
             ca06(n) = acof(n,5,mt) - acof(n,7,mt)
             ca07(n) = acof(n,4,mt) + acof(n,8,mt)
             ca08(n) = acof(n,6,mt)
             !
             !        f-states: a cubic set (cub)
             !
             ca09(n) = ( acof(n,9,mt)  - acof(n,15,mt) )*c5 -&
                  ( acof(n,11,mt) - acof(n,13,mt) )*c3
             ca10(n) = ( acof(n,9,mt)  + acof(n,15,mt) )*c5 +&
                  ( acof(n,11,mt) + acof(n,13,mt) )*c3
             ca11(n) =   acof(n,12,mt)
             ca12(n) = ( acof(n,9,mt)  + acof(n,15,mt) )*c3 -&
                  ( acof(n,11,mt) + acof(n,13,mt) )*c5
             ca13(n) =   acof(n,10,mt) + acof(n,14,mt)
             ca14(n) = ( acof(n,9,mt)  - acof(n,15,mt) )*c3 +&
                  ( acof(n,11,mt) - acof(n,13,mt) )*c5
             ca15(n) =   acof(n,10,mt) - acof(n,14,mt)

             !
             !bcof
             !       s-states
             cb00(n) =  bcof(n,0,mt)
             !       p-states
             cb01(n) =  bcof(n,1,mt) - bcof(n,3,mt)
             cb02(n) =  bcof(n,1,mt) + bcof(n,3,mt)
             cb03(n) =  bcof(n,2,mt)
             !        d-states
             cb04(n) =  bcof(n,4,mt) - bcof(n,8,mt)
             cb05(n) =  bcof(n,5,mt) + bcof(n,7,mt)
             cb06(n) =  bcof(n,5,mt) - bcof(n,7,mt)
             cb07(n) =  bcof(n,4,mt) + bcof(n,8,mt)
             cb08(n) =  bcof(n,6,mt)
             !
             !       f-states: a cubic set (cub)
             !
             cb09(n) = ( bcof(n,9,mt)  - bcof(n,15,mt) )*c5 -&
                  ( bcof(n,11,mt) - bcof(n,13,mt) )*c3
             cb10(n) = ( bcof(n,9,mt)  + bcof(n,15,mt) )*c5 +&
                  ( bcof(n,11,mt) + bcof(n,13,mt) )*c3
             cb11(n) =   bcof(n,12,mt)
             cb12(n) = ( bcof(n,9,mt)  + bcof(n,15,mt) )*c3 -&
                  ( bcof(n,11,mt) + bcof(n,13,mt) )*c5
             cb13(n) =   bcof(n,10,mt) + bcof(n,14,mt)
             cb14(n) = ( bcof(n,9,mt)  - bcof(n,15,mt) )*c3 +&
                  ( bcof(n,11,mt) - bcof(n,13,mt) )*c5
             cb15(n) =   bcof(n,10,mt) - bcof(n,14,mt)
          ENDDO !n

          DO n=1,SIZE(acof,1)
             DO m=1,SIZE(acof,1)
                !  s
                mmn(1,orbi,m,n)  =   ca00(n)*CONJG(ca00(m)) &
                     + cb00(n)*CONJG(cb00(m))*ddn0
                !  p
                mmn(2,orbi,m,n)  = ( ca01(n)*CONJG(ca01(m)) &
                     + cb01(n)*CONJG(cb01(m))*ddn1 )*h
                mmn(3,orbi,m,n)  = ( ca02(n)*CONJG(ca02(m)) &
                     + cb02(n)*CONJG(cb02(m))*ddn1 )*h
                mmn(4,orbi,m,n)  =   ca03(n)*CONJG(ca03(m)) &
                     + cb03(n)*CONJG(cb03(m))*ddn1
                !  d
                mmn(5,orbi,m,n)  = ( ca04(n)*CONJG(ca04(m)) &
                     + cb04(n)*CONJG(cb04(m))*ddn2 )*h
                mmn(6,orbi,m,n)  = ( ca05(n)*CONJG(ca05(m)) &
                     + cb05(n)*CONJG(cb05(m))*ddn2 )*h
                mmn(7,orbi,m,n)  = ( ca06(n)*CONJG(ca06(m)) &
                     + cb06(n)*CONJG(cb06(m))*ddn2 )*h
                mmn(8,orbi,m,n)  = ( ca07(n)*CONJG(ca07(m)) &
                     + cb07(n)*CONJG(cb07(m))*ddn2 )*h
                mmn(9,orbi,m,n)  =   ca08(n)*CONJG(ca08(m)) &
                     + cb08(n)*CONJG(cb08(m))*ddn2
                !  f: a cubic set
                IF(l_oc_f)THEN
                   mmn(10,orbi,m,n) = ( ca09(n)*CONJG(ca09(m)) &
                        + cb09(n)*CONJG(cb09(m))*ddn3 )*g
                   mmn(11,orbi,m,n) = ( ca10(n)*CONJG(ca10(m)) &
                        + cb10(n)*CONJG(cb10(m))*ddn3 )*g
                   mmn(12,orbi,m,n) =   ca11(n)*CONJG(ca11(m)) &
                        + cb11(n)*CONJG(cb11(m))*ddn3
                   mmn(13,orbi,m,n) = ( ca12(n)*CONJG(ca12(m)) &
                        + cb12(n)*CONJG(cb12(m))*ddn3 )*g
                   mmn(14,orbi,m,n) = ( ca13(n)*CONJG(ca13(m)) &
                        + cb13(n)*CONJG(cb13(m))*ddn3 )*h
                   mmn(15,orbi,m,n) = ( ca14(n)*CONJG(ca14(m)) &
                        + cb14(n)*CONJG(cb14(m))*ddn3 )*g
                   mmn(16,orbi,m,n) = ( ca15(n)*CONJG(ca15(m)) &
                        + cb15(n)*CONJG(cb15(m))*ddn3 )*h
                ENDIF !l_oc_f
             ENDDO !m
          ENDDO !n

          !--------------------------------------------------------------------
          !ccof   ( contributions from local orbitals )
          !
          lo_loop:DO lo = 1,atoms%nlo(ityp)
             l = atoms%llo(lo,ityp)
             !lo-s
             IF ( l.EQ.0 )  THEN
                DO n=1,SIZE(acof,1)
                   cc00(n) = ccof(0,n,lo,mt)
                ENDDO
                !
                DO n=1,SIZE(acof,1)
                   DO m=1,SIZE(acof,1)
                      mmn(1,orbi,m,n)=mmn(1,orbi,m,n)+&
                           ( ca00(n)*CONJG(cc00(m)) + &
                           cc00(n)*CONJG(ca00(m)) )*usdus%uulon(lo,ityp,jspin) + &
                           ( cb00(n)*CONJG(cc00(m)) + &
                           cc00(n)*CONJG(cb00(m)) )*usdus%dulon(lo,ityp,jspin) +&
                           cc00(n)*CONJG(cc00(m))*usdus%uloulopn(lo,lo,ityp,jspin)
                   ENDDO !m
                ENDDO !n
                CYCLE lo_loop
             ENDIF
             !lo-p
             IF ( l.EQ.1 )  THEN
                DO n=1,SIZE(acof,1)
                   cc01(n) = ccof(-1,n,lo,mt) - ccof(1,n,lo,mt)
                   cc02(n) = ccof(-1,n,lo,mt) + ccof(1,n,lo,mt)
                   cc03(n) = ccof( 0,n,lo,mt)
                ENDDO
                DO n=1,SIZE(acof,1)
                   DO m=1,SIZE(acof,1)
                      mmn(2,orbi,m,n) = mmn(2,orbi,m,n)  +&
                           (   ( ca01(n)*CONJG(cc01(m)) + &
                           cc01(n)*CONJG(ca01(m)) )*usdus%uulon(lo,ityp,jspin) +&
                           ( cb01(n)*CONJG(cc01(m)) + &
                           cc01(n)*CONJG(cb01(m)) )*usdus%dulon(lo,ityp,jspin) +&
                           cc01(n)*CONJG(cc01(m))*usdus%uloulopn(lo,lo,ityp,jspin) )*h

                      mmn(3,orbi,m,n) = mmn(3,orbi,m,n)  +&
                           (   ( ca02(n)*CONJG(cc02(m)) + &
                           cc02(n)*CONJG(ca02(m)) )*usdus%uulon(lo,ityp,jspin) +&
                           ( cb02(n)*CONJG(cc02(m)) + &
                           cc02(n)*CONJG(cb02(m)) )*usdus%dulon(lo,ityp,jspin) +&
                           cc02(n)*CONJG(cc02(m))*usdus%uloulopn(lo,lo,ityp,jspin) )*h

                      mmn(4,orbi,m,n) = mmn(4,orbi,m,n)  +&
                           ( ca03(n)*CONJG(cc03(m)) + &
                           cc03(n)*CONJG(ca03(m)) )*usdus%uulon(lo,ityp,jspin) +&
                           ( cb03(n)*CONJG(cc03(m)) + &
                           cc03(n)*CONJG(cb03(m)) )*usdus%dulon(lo,ityp,jspin) +&
                           cc03(n)*CONJG(cc03(m))*usdus%uloulopn(lo,lo,ityp,jspin)

                   ENDDO !m
                ENDDO !n
                CYCLE lo_loop
             ENDIF
             !lo-d
             IF ( l.EQ.2 )  THEN
                DO n=1,SIZE(acof,1)
                   cc04(n) = ccof(-2,n,lo,mt) - ccof(2,n,lo,mt)
                   cc05(n) = ccof(-1,n,lo,mt) + ccof(1,n,lo,mt)
                   cc06(n) = ccof(-1,n,lo,mt) - ccof(1,n,lo,mt)
                   cc07(n) = ccof(-2,n,lo,mt) + ccof(2,n,lo,mt)
                   cc08(n) = ccof( 0,n,lo,mt)
                ENDDO !n
                !
                DO n=1,SIZE(acof,1)
                   DO m=1,SIZE(acof,1)
                      mmn(5,orbi,m,n) = mmn(5,orbi,m,n)  +&
                           (( ca04(n)*CONJG(cc04(m)) + &
                           cc04(n)*CONJG(ca04(m)) )*usdus%uulon(lo,ityp,jspin) +&
                           ( cb04(n)*CONJG(cc04(m)) + &
                           cc04(n)*CONJG(cb04(m)) )*usdus%dulon(lo,ityp,jspin) +&
                           cc04(n)*CONJG(cc04(m))*usdus%uloulopn(lo,lo,ityp,jspin))*h
                      mmn(6,orbi,m,n) = mmn(6,orbi,m,n)  +&
                           (( ca05(n)*CONJG(cc05(m)) + &
                           cc05(n)*CONJG(ca05(m)) )*usdus%uulon(lo,ityp,jspin) +&
                           ( cb05(n)*CONJG(cc05(m)) + &
                           cc05(n)*CONJG(cb05(m)) )*usdus%dulon(lo,ityp,jspin) +&
                           cc05(n)*CONJG(cc05(m))*usdus%uloulopn(lo,lo,ityp,jspin) )*h
                      mmn(7,orbi,m,n) = mmn(7,orbi,m,n)  +&
                           (( ca06(n)*CONJG(cc06(m)) + &
                           cc06(n)*CONJG(ca06(m)) )*usdus%uulon(lo,ityp,jspin) +&
                           ( cb06(n)*CONJG(cc06(m)) + &
                           cc06(n)*CONJG(cb06(m)) )*usdus%dulon(lo,ityp,jspin) +&
                           cc06(n)*CONJG(cc06(m))*usdus%uloulopn(lo,lo,ityp,jspin) )*h
                      mmn(8,orbi,m,n) = mmn(8,orbi,m,n)  +&
                           (( ca07(n)*CONJG(cc07(m)) + &
                           cc07(n)*CONJG(ca07(m)) )*usdus%uulon(lo,ityp,jspin) +&
                           ( cb07(n)*CONJG(cc07(m)) + &
                           cc07(n)*CONJG(cb07(m)) )*usdus%dulon(lo,ityp,jspin) +&
                           cc07(n)*CONJG(cc07(m))*usdus%uloulopn(lo,lo,ityp,jspin) )*h
                      mmn(9,orbi,m,n) = mmn(9,orbi,m,n)  +&
                           ( ca08(n)*CONJG(cc08(m)) + &
                           cc08(n)*CONJG(ca08(m)) )*usdus%uulon(lo,ityp,jspin) +&
                           ( cb08(n)*CONJG(cc08(m)) + &
                           cc08(n)*CONJG(cb08(m)) )*usdus%dulon(lo,ityp,jspin) +&
                           cc08(n)*CONJG(cc08(m))*usdus%uloulopn(lo,lo,ityp,jspin)
                   ENDDO !m
                ENDDO !n
                CYCLE lo_loop
             ENDIF
             !lo-f
             IF ( l.EQ.3 .AND. l_oc_f)  THEN
                !
                !  a cubic set (cub)
                !
                DO n=1,SIZE(acof,1)
                   cc09(n) = ( ccof(-3,n,lo,mt) - ccof(3,n,lo,mt) )*c5 -&
                        ( ccof(-1,n,lo,mt) - ccof(1,n,lo,mt) )*c3
                   cc10(n) = ( ccof(-3,n,lo,mt) + ccof(3,n,lo,mt) )*c5 +&
                        ( ccof(-1,n,lo,mt) + ccof(1,n,lo,mt) )*c3
                   cc11(n) =   ccof( 0,n,lo,mt)
                   cc12(n) = ( ccof(-3,n,lo,mt) + ccof(3,n,lo,mt) )*c3 -&
                        ( ccof(-1,n,lo,mt) + ccof(1,n,lo,mt) )*c5
                   cc13(n) =   ccof(-2,n,lo,mt) + ccof(2,n,lo,mt)
                   cc14(n) = ( ccof(-3,n,lo,mt) - ccof(3,n,lo,mt) )*c3 +&
                        ( ccof(-1,n,lo,mt) - ccof(1,n,lo,mt) )*c5
                   cc15(n) =   ccof(-2,n,lo,mt) - ccof(2,n,lo,mt)
                ENDDO !n
                DO n=1,SIZE(acof,1)
                   DO m=1,SIZE(acof,1)
                      mmn(10,orbi,m,n) = mmn(10,orbi,m,n)  +&
                           (( ca09(n)*CONJG(cc09(m)) + &
                           cc09(n)*CONJG(ca09(m)) )*usdus%uulon(lo,ityp,jspin) +&
                           ( cb09(n)*CONJG(cc09(m)) + &
                           cc09(n)*CONJG(cb09(m)) )*usdus%dulon(lo,ityp,jspin) +&
                           cc09(n)*CONJG(cc09(m))*usdus%uloulopn(lo,lo,ityp,jspin) )*g
                      mmn(11,orbi,m,n) = mmn(11,orbi,m,n)  +&
                           (( ca10(n)*CONJG(cc10(m)) + &
                           cc10(n)*CONJG(ca10(m)) )*usdus%uulon(lo,ityp,jspin) +&
                           ( cb10(n)*CONJG(cc10(m)) + &
                           cc10(n)*CONJG(cb10(m)) )*usdus%dulon(lo,ityp,jspin) +&
                           cc10(n)*CONJG(cc10(m))*usdus%uloulopn(lo,lo,ityp,jspin) )*g
                      mmn(12,orbi,m,n) = mmn(12,orbi,m,n)  +&
                           ( ca11(n)*CONJG(cc11(m)) + &
                           cc11(n)*CONJG(ca11(m)) )*usdus%uulon(lo,ityp,jspin) +&
                           ( cb11(n)*CONJG(cc11(m)) + &
                           cc11(n)*CONJG(cb11(m)) )*usdus%dulon(lo,ityp,jspin) +&
                           cc11(n)*CONJG(cc11(m))*usdus%uloulopn(lo,lo,ityp,jspin)
                      mmn(13,orbi,m,n) = mmn(13,orbi,m,n)  +&
                           (( ca12(n)*CONJG(cc12(m)) + &
                           cc12(n)*CONJG(ca12(m)) )*usdus%uulon(lo,ityp,jspin) +&
                           ( cb12(n)*CONJG(cc12(m)) + &
                           cc12(n)*CONJG(cb12(m)) )*usdus%dulon(lo,ityp,jspin) +&
                           cc12(n)*CONJG(cc12(m))*usdus%uloulopn(lo,lo,ityp,jspin) )*g
                      mmn(14,orbi,m,n) = mmn(14,orbi,m,n)  +&
                           (( ca13(n)*CONJG(cc13(m)) + &
                           cc13(n)*CONJG(ca13(m)) )*usdus%uulon(lo,ityp,jspin) +&
                           ( cb13(n)*CONJG(cc13(m)) + &
                           cc13(n)*CONJG(cb13(m)) )*usdus%dulon(lo,ityp,jspin) +&
                           cc13(n)*CONJG(cc13(m))*usdus%uloulopn(lo,lo,ityp,jspin) )*h
                      mmn(15,orbi,m,n) = mmn(15,orbi,m,n)  +&
                           (( ca14(n)*CONJG(cc14(m)) + &
                           cc14(n)*CONJG(ca14(m)) )*usdus%uulon(lo,ityp,jspin) +&
                           ( cb14(n)*CONJG(cc14(m)) + &
                           cc14(n)*CONJG(cb14(m)) )*usdus%dulon(lo,ityp,jspin) +&
                           cc14(n)*CONJG(cc14(m))*usdus%uloulopn(lo,lo,ityp,jspin) )*g
                      mmn(16,orbi,m,n) = mmn(16,orbi,m,n)  +&
                           (( ca15(n)*CONJG(cc15(m)) + &
                           cc15(n)*CONJG(ca15(m)) )*usdus%uulon(lo,ityp,jspin) +&
                           ( cb15(n)*CONJG(cc15(m)) + &
                           cc15(n)*CONJG(cb15(m)) )*usdus%dulon(lo,ityp,jspin) +&
                           cc15(n)*CONJG(cc15(m))*usdus%uloulopn(lo,lo,ityp,jspin) )*h
                   ENDDO !m
                ENDDO !n
             ENDIF
          ENDDO lo_loop

       ENDDO !rep

    ENDDO     ! types (ityp)

  END SUBROUTINE wann_orbcomp
END MODULE m_wann_orbcomp

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


subroutine sub_measure(vars)
  use mod_variables,only : imax,jmax,imx1,jmx1,n_prb,n_mes,variances,  &
                           xcen,ycen,ih,ihx,ihy,iskip,jskip,mbot,mtop, &
                           mlft,mrht
  implicit none
!
  type(variances) :: vars
  integer :: ic,jc,indx
  integer :: iht(imax,jmax)
!
!
  iht(:,:) = 0
!
  do jc=1,jmax,iskip
  do ic=1,imax,jskip
    if(mlft<xcen(ic) .and. xcen(ic)<mrht .and. mbot<ycen(jc) .and. ycen(jc)<mtop) then
      iht(ic,jc) = 1
    else
      iht(ic,jc) = 0
    endif
  enddo
  enddo
!
  allocate(ihx(imx1,jmax),ihy(imax,jmx1))
  ihx(:,:) = 0
  ihy(:,:) = 0
!
  do jc=1,jmax
  do ic=1,imax
    if(iht(ic,jc)==1) then
      ihx(ic  ,jc) = 1
      ihx(ic+1,jc) = 1
      ihy(ic,jc  ) = 1
      ihy(ic,jc+1) = 1
    endif
  enddo
  enddo
  n_mes = 0
  do jc=1,jmax
  do ic=1,imx1
    if(ihx(ic,jc)==1) then
      n_mes = n_mes+1
    endif
  enddo
  enddo
  do jc=1,jmx1
  do ic=1,imax
    if(ihy(ic,jc)==1) then
      n_mes = n_mes+1
    endif
  enddo
  enddo
!
!
  allocate(ih(n_mes))
  ih(:) = 0
  n_mes = 0
  indx  = 0
  do jc=1,jmax
  do ic=1,imx1
    indx = indx+1
    if(ihx(ic,jc)==1) then
      n_mes = n_mes+1
      ih(n_mes) = indx
    endif
  enddo
  enddo
  do jc=1,jmx1
  do ic=1,imax
    indx = indx+1
    if(ihy(ic,jc)==1) then
      n_mes = n_mes+1
      ih(n_mes) = indx
    endif
  enddo
  enddo
  write(*,*)
  write(*,*) "Number of measurement points:",n_mes
!
  return
end subroutine sub_measure


text:
  subroutine genqsat(temp,qsat,lat2,lon2,ps,icesat,slg,a5,b5)
!   input argument list:
!     t,ps     - temperature, psfc on gaussian grid
!     lat1     - number of gaussian lats in the sub-domain array
!     lon1     - number of gaussian longitudes in the sub-domain array
!     nsig     - number of sigma levels
!     sigl      - sigma values at mid-point of model levels
!     icesat   - logical flag:  T=include ice and ice-water effects,
!                depending on t, in qsat and esat calcuations.
!                otherwise, compute esat and qsat with respect to water surface
!                depending on t, in qsat and esat calcuations.
!                otherwise, compute esat and qsat with respect to water surface
!
!   output argument list:
!     qsat     - specific humidity (input), saturation specific humidity (output)
  use type_kinds, only: fp_kind
  use variables,only: hybrid,nsig,rd,zero,psat,fv,one,rv,cliq,cvap,hvap,csol,&
      hfus,ttp
  implicit none

  real(fp_kind) dldt,xa,xb,hsub,dldti,xai,xbi,tmix
  parameter(dldt=cvap-cliq,xa=-dldt/rv,xb=xa+hvap/(rv*ttp))
  parameter(hsub=hvap+hfus)
  parameter(dldti=cvap-csol,xai=-dldti/rv,xbi=xai+hsub/(rv*ttp))
  parameter (tmix=ttp-20.0_fp_kind)
  real(fp_kind) eps,omeps
  parameter (eps=rd/rv)
  parameter (omeps = 1._fp_kind-eps)
!
  logical icesat
  integer k,j,i,ntsig
  integer lat2,lon2
  real(fp_kind) pw,q,tdry,tr,es,qs,esi,esw
  real(fp_kind) w,onep3,esmax
  real(fp_kind),dimension(lat2,lon2,nsig):: qsat
! real(fp_kind),dimension(lat2,lon2,nsig):: esat
  real(fp_kind),dimension(lat2,lon2,nsig):: temp
  real(fp_kind),dimension(lat2,lon2):: ps
  real(fp_kind),dimension(nsig):: slg
  real(fp_kind),dimension(nsig+1):: a5,b5
!
  onep3 = 1.e3_fp_kind

  if (icesat) then
     do k = 1,nsig
        do j = 1,lon2
           do i = 1,lat2

             if (hybrid) then
               pw = 0.5*(a5(k)+a5(k+1)+(b5(k)+b5(k+1))* &
                          ps(i,j))
             else
               pw = ps(i,j)*slg(k)
             end if
              pw  = onep3*pw
! maximum vapor pressure 5% of atmospheric pressure
              esmax=0.05*pw

              q  = qsat(i,j,k)
              if (q.lt.zero) q=zero

              tdry = temp(i,j,k)/(one+fv*q)
              tr = ttp/tdry
              if (tdry >= ttp) then
                 es = psat * (tr**xa) * exp(xb*(one-tr))
              elseif (tdry < tmix) then
                 es = psat * (tr**xai) * exp(xbi*(one-tr))
              else
                 w  = (tdry - tmix) / (ttp - tmix)
                 es =  w * psat * (tr**xa) * exp(xb*(one-tr)) &
                      + (one-w) * psat * (tr**xai) * exp(xbi*(one-tr))
              endif

              es = min(es,esmax)
              qs = eps * es / (pw - omeps * es)

              if (qs.lt.qsat(i,j,k)) then
                 tdry = temp(i,j,k)/(one+fv*qs)
                 tr = ttp/tdry
                 if (tdry >= ttp) then
                    es = psat * (tr**xa) * exp(xb*(one-tr))
                 elseif (tdry < tmix) then
                    es = psat * (tr**xai) * exp(xbi*(one-tr))
                 else
                    w  = (tdry - tmix) / (ttp - tmix)
                    es =  w * psat * (tr**xa) * exp(xb*(one-tr)) &
                         + (one-w) * psat * (tr**xai) * exp(xbi*(one-tr))
                 endif

                 es = min(es,esmax)
                 qs = eps * es / (pw - omeps * es)

                 tdry = temp(i,j,k)/(one+fv*qs)
                 tr = ttp/tdry
                 if (tdry >= ttp) then
                    es = psat * (tr**xa) * exp(xb*(one-tr))
                 elseif (tdry < tmix) then
                    es = psat * (tr**xai) * exp(xbi*(one-tr))
                 else
                    w  = (tdry - tmix) / (ttp - tmix)
                    es =  w * psat * (tr**xa) * exp(xb*(one-tr)) &
                         + (one-w) * psat * (tr**xai) * exp(xbi*(one-tr))
                 endif
                 es = min(es,esmax)
                 qs = eps * es / (pw - omeps * es)

              end if

              qsat(i,j,k) = qs

           end do
        end do
     end do

!
!     Compute saturation values with respect to water surface
  else
     do k = 1,nsig
        do j = 1,lon2
           do i = 1,lat2

             if (hybrid) then
               pw = 0.5*(a5(k)+a5(k+1)+(b5(k)+b5(k+1))* &
                          ps(i,j))
             else
               pw = ps(i,j)*slg(k)
             end if

              pw  = onep3*pw
! maximum vapor pressure 5% of atmospheric pressure
              esmax=0.05*pw

              q  = qsat(i,j,k)
              if (q.lt.zero) q=zero

              tdry = temp(i,j,k)/(one+fv*q)
              tr = ttp/tdry
              es = psat * (tr**xa) * exp(xb*(one-tr))
              es = min(es,esmax)
              qs = eps * es / (pw - omeps * es)

              if (qs.lt.qsat(i,j,k)) then
                 tdry = temp(i,j,k)/(one+fv*qs)
                 tr = ttp/tdry
                 es = psat * (tr**xa) * exp(xb*(one-tr))
                 es = min(es,esmax)
                 qs = eps * es / (pw - omeps * es)
                 tdry = temp(i,j,k)/(one+fv*qs)
                 tr = ttp/tdry
                 es = psat * (tr**xa) * exp(xb*(one-tr))
                 es = min(es,esmax)
                 qs = eps * es / (pw - omeps * es)

              end if
              qsat(i,j,k) = qs

           end do
        end do
     end do

  endif

  return
  end subroutine genqsat

max_stars_repo_name: GEOS-ESM/GEOSadas

max_issues_repo_name: GEOS-ESM/GEOSadas

max_forks_repo_name: GEOS-ESM/GEOSadas


PROGRAM LA_SSBGVX_ET_EXAMPLE
!
!  -- LAPACK95 interface driver routine (version 3.0) --
!     UNI-C, Denmark; Univ. of Tennessee, USA; NAG Ltd., UK
!     September, 2000
!
!  .. USE STATEMENTS
   USE LA_PRECISION, ONLY: WP => SP
   USE F95_LAPACK, ONLY: LA_SBGVX
!  .. IMPLICIT STATEMENT ..
   IMPLICIT NONE
!  .. PARAMETERS ..
   CHARACTER(LEN=*), PARAMETER :: FMT = '(8(1X,F10.3))'
   INTEGER, PARAMETER :: NIN=5, NOUT=6
!  .. LOCAL SCALARS ..
   CHARACTER(LEN=1) :: UPLO
   INTEGER :: I, J, INFO, N, KA, KB, M, IL, IU
   REAL(WP) ::  VL, VU
!  .. LOCAL ARRAYS ..
   INTEGER, ALLOCATABLE :: IFAIL(:)
   REAL(WP), ALLOCATABLE :: AA(:,:), BB(:,:), W(:)
   REAL(WP), ALLOCATABLE :: A(:,:), B(:,:), Z(:,:), DUMMY(:,:), Q(:,:)
!  .. EXECUTABLE STATEMENTS ..
   WRITE (NOUT,*) 'SSBGVX ET_Example Program Results.'
   READ ( NIN, * )   ! SKIP HEADING IN DATA FILE
   READ ( NIN, * ) N, KA, KB
   PRINT *, 'N = ', N, ' KA = ', KA, ' KB = ', KB
   ALLOCATE ( A(KA+1,N), AA(KA+1,N), B(KB+1,N), BB(KB+1,N), W(N), Z(N,N) )
   ALLOCATE (IFAIL(N), Q(N,N))
!
   VL = -10.0; VU=10.0; IL=1; IU=N; UPLO = 'U'
   AA = HUGE(1.0_WP); BB = HUGE(1.0_WP)
   DO I = 1, KA+1; READ (NIN, *) (AA(I, J), J = KA-I+2, N); ENDDO
   DO I = 1, KB+1; READ (NIN, *) (BB(I, J), J = KB-I+2, N); ENDDO
   A=AA; B=BB
   WRITE(NOUT,*) 'The matrix A:'
   DO I = 1, KA+1; WRITE (NOUT,*) 'I = ', I; WRITE (NOUT,FMT) AA(I,:); ENDDO
   WRITE(NOUT,*) 'The matrix B:'
   DO I = 1, KB+1; WRITE (NOUT,*) 'I = ', I; WRITE (NOUT,FMT) BB(I,:); ENDDO
!
   WRITE ( NOUT, * )'---------------------------------------------------------'
   WRITE ( NOUT, * )
   WRITE ( NOUT, * )'Details of LA_SSBGVX LAPACK Subroutine Results.'
   WRITE ( NOUT, * )
!
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_SBGVX( A, B, W, ''U'', Z, INFO )'
   A=AA; B=BB
   CALL LA_SBGVX( A, B, W, UPLO, Z, INFO = INFO )
   WRITE(NOUT,*) 'INFO = ', INFO, ' EIGENVALUES:'
   WRITE(NOUT,FMT) W
   WRITE(NOUT,*) 'EIGENVECTORS:'
   DO I = 1, N; WRITE(NOUT,*) 'I = ', I; WRITE (NOUT,FMT) Z(:,I); END DO
!
     UPLO = 'U'
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_SBGVX( A, B, W, ''U'', Z )'
   A=AA; B=BB
   CALL LA_SBGVX( A, B, W, UPLO, Z )
   WRITE(NOUT,*) ' EIGENVALUES:'
   WRITE(NOUT,FMT) W
   WRITE(NOUT,*) 'EIGENVECTORS:'
   DO I = 1, N; WRITE(NOUT,*) 'I = ', I; WRITE (NOUT,FMT) Z(:,I); END DO
!
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_SBGVX( A, B, W, Z=Z )'
   A=AA; B=BB
   CALL LA_SBGVX( A, B, W, Z=Z )
   WRITE(NOUT,*) ' EIGENVALUES:'
   WRITE(NOUT,FMT) W
   WRITE(NOUT,*) 'EIGENVECTORS:'
   DO I = 1, N; WRITE(NOUT,*) 'I = ', I; WRITE (NOUT,FMT) Z(:,I); END DO
!
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_SBGVX( A, B, W, ''L'', Z, INFO=INFO )'
   A=AA; B=BB; W = HUGE(1.0_WP); Z = HUGE(1.0_WP)
   CALL LA_SBGVX( A, B, W, UPLO, Z, INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO, ' EIGENVALUES:'
   WRITE(NOUT,FMT) W
   WRITE(NOUT,*) 'EIGENVECTORS:'
   DO I = 1, N; WRITE(NOUT,*) 'I = ', I; WRITE (NOUT,FMT) Z(:,I); END DO
!
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_SBGVX( A, B, W )'
   A=AA; B=BB; Z = HUGE(1.0_WP)
   CALL LA_SBGVX( A, B, W )
   WRITE(NOUT,*) ' EIGENVALUES:'
   WRITE(NOUT,FMT) W
! STARTING THE ERROR TESTS:
! ERROR 1
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_SBGVX( DUMMY, B, W, INFO=INFO )'
   A=AA; B=BB; Z = HUGE(1.0_WP)
   CALL LA_SBGVX( DUMMY, B, W, INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
! ERROR 2
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_SBGVX( A, DUMMY, W, INFO=INFO )'
   A=AA; B=BB; Z = HUGE(1.0_WP)
   CALL LA_SBGVX( A, DUMMY, W, INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
! ERROR 3
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_SBGVX( A, B, W(1:N-1), INFO=INFO )'
   A=AA; B=BB; Z = HUGE(1.0_WP)
   CALL LA_SBGVX( A, B, W(1:N-1), INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
! ERROR 4
   UPLO = '9'
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_SBGVX( A, B, W, UPLO=''9'', INFO=INFO )'
   A=AA; B=BB; Z = HUGE(1.0_WP)
   CALL LA_SBGVX( A, B, W, UPLO, INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
! ERROR 5
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_SBGVX( A, B, W, UPLO, Z(:, 1:N-1), INFO=INFO )'
   A=AA; B=BB; Z = HUGE(1.0_WP)
   CALL LA_SBGVX( A, B, W, UPLO, Z(:, 1:N-1), INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
! ERROR 5
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_SBGVX( A, B, W, UPLO, Z(1:N-1,:), INFO=INFO )'
   A=AA; B=BB; Z = HUGE(1.0_WP)
   CALL LA_SBGVX( A, B, W, UPLO, Z(1:N-1,:), INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
! ERROR 6
   VL = 10.0; VU = -10.0
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_SBGVX( A, B, W, UPLO, Z, VL, VU, M, IFAIL, Q, INFO=INFO )'
   A=AA; B=BB; Z = HUGE(1.0_WP)
   CALL LA_SBGVX( A, B, W, UPLO, Z,  VL, VU, M=M, IFAIL=IFAIL, Q=Q, INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
! ERROR 7
   IL = 1; IU = N; VL=-10.0; VU=10.0
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_SBGVX( A, B, W, UPLO, Z, VL, VU, IL, IU, M, IFAIL, Q, INFO=INFO )'
   A=AA; B=BB; Z = HUGE(1.0_WP)
   CALL LA_SBGVX( A, B, W, UPLO, Z, VL, VU, IL, IU, M, IFAIL, Q, INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
   IU=N; IL=1
! ERROR 8
   IL = -1; IU = N+1
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_SBGVX( A, B, W, UPLO, Z, IL=IL, IU=IU, M=M, IFAIL=IFAIL, Q=Q, INFO=INFO )'
   A=AA; B=BB; Z = HUGE(1.0_WP)
   CALL LA_SBGVX( A, B, W, UPLO, Z, IL=IL, IU=IU, M=M, IFAIL=IFAIL, Q=Q, INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
   IU=N; IL=1
! ERROR 9
   IU = N+1
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_SBGVX( A, B, W, UPLO, Z, IL=IL, IU=IU, M=M, IFAIL=IFAIL, Q=Q, INFO=INFO )'
   A=AA; B=BB; Z = HUGE(1.0_WP)
   CALL LA_SBGVX( A, B, W, UPLO, Z, IL=IL, IU=IU, M=M, IFAIL=IFAIL, Q=Q, INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
   IU=N
! ERROR 10
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_SBGVX( A, B, W,  UPLO, Z, VL, VU, M=M, IFAIL=IFAIL(1:N-1), Q=Q, INFO=INFO )'
   A=AA; B=BB; Z = HUGE(1.0_WP)
   CALL LA_SBGVX( A, B, W, UPLO, Z, VL, VU, M=M, IFAIL=IFAIL(1:N-1), Q=Q, INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
! ERROR 10
   UPLO= 'L'
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_SBGVX( A, B, W,  UPLO, VL=VL, VU=VU, M=M, IFAIL=IFAIL, Q=Q, INFO=INFO )'
   A=AA; B=BB; Z = HUGE(1.0_WP)
   CALL LA_SBGVX( A, B, W, UPLO, VL=VL, VU=VU, M=M, IFAIL=IFAIL, Q=Q, INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
! ERROR 11
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_SBGVX( A, B, W, UPLO, Z, VL, VU, M=M, IFAIL=IFAIL, Q=Q(1:N-1,:), INFO=INFO )'
   A=AA; B=BB; Z = HUGE(1.0_WP)
   CALL LA_SBGVX( A, B, W, UPLO, Z,  VL, VU, M=M, IFAIL=IFAIL, Q=Q(1:N-1,:), INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
! ERROR 11
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_SBGVX( A, B, W, UPLO, Z, VL, VU, M=M, IFAIL=IFAIL, Q=Q(:,1:N-1), INFO=INFO )'
   A=AA; B=BB; Z = HUGE(1.0_WP)
   CALL LA_SBGVX( A, B, W, UPLO, Z,  VL, VU, M=M, IFAIL=IFAIL, Q=Q(:,1:N-1), INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
! ERROR 11
   WRITE(NOUT,*)
   WRITE(NOUT,*) 'CALL LA_SBGVX( A, B, W, UPLO, VL=VL, VU=VU, M=M, IFAIL=IFAIL, Q=Q, INFO=INFO )'
   A=AA; B=BB; Z = HUGE(1.0_WP)
   CALL LA_SBGVX( A, B, W, UPLO, VL=VL, VU=VU, M=M, IFAIL=IFAIL, Q=Q, INFO=INFO )
   WRITE(NOUT,*) 'INFO = ', INFO
!
END PROGRAM LA_SSBGVX_ET_EXAMPLE


text: #redirect California Aggie Marching BandUh

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
#if defined(CRAY)
#define ERFC derfc
#else
#define ERFC derfc
#endif
      subroutine find_Ewald_eta_screened(rb, errlim, nat,
     $     atp, eps, eta)
c=================================================================
c     Calculates optimal "eta" for Ewald summations.
c
c   INPUT:
c     rb - real-space basis
c     errlim - requested accuracy of summations
c     nat - number of atoms
c     atp - atomic positions
c     eps - high-frequency dielectric tensor
c
c   OUTPUT:
c     eta - parameter for the Ewald summations
c
c   DATE: Sun Jan  1 12:08:58 PST 2012
c=================================================================
      implicit none

      integer MAXIT
      parameter( MAXIT = 100 )

      integer nat
      real*8 rb(3,3), eta, errlim, atp(3,nat), eps(3,3)

      integer nq1, nq2, nq3, nr1, nr2, nr3, m, n, iter, j
      real*8 dtau(3), rsq, taumax, ekmax, pi, tpi, tpi2,
     $     qb(3,3), x, y, FACT, etaprev

c     LAPACK related variables
      integer lwork
      parameter (lwork=30)
      integer info
      real*8 teps(3,3), evr(3), evi(3), vleft(3,3), vright(3,3),
     $     work(lwork)

c     Complementary error function.
      real*8, external                          :: ERFC

!      print *
      if(errlim > 1d0) then
         print *, 'errlim in find_Ewald_eta is too large = ', errlim
         stop
      end if

      pi = 4d0*atan(1d0)
      tpi = 2d0*pi
      tpi2 = tpi**2

      teps = eps
      call DGEEV('N','V',3,teps,3,evr,evi,vleft,3,vright,3,
     $     work,lwork,info)
      if( info /= 0 ) print *, 'DGEEV error: info=', info

#ifdef DEBUG
      print '(2x,a30)', 'Eigenvalues and eigenvectors of ',
     $     'dielectric matrix:'
      do j = 1, 3
         print '(2f10.4,5x,3f10.4)', evr(j), evi(j), vright(:,j)
      end do
#endif

      call lattc(rb,qb)
      taumax = 0d0
      do m=1,nat
         do n=1,m
            dtau = atp(:,m) - atp(:,n)
            rsq = dot_product(dtau, dtau)
            taumax = max(taumax, rsq)
         end do
      end do

      if( abs(eta) < 1d-5 ) eta = 1d0
      FACT = 2d0
      etaprev = eta

      do iter = 1, MAXIT
         x = 2d0 * eta * sqrt(-log(errlim))
         ekmax = MINVAL(1d0/evr) * x**2 / 2d0
         call gbox(ekmax, 1d0, qb, nq1, nq2, nq3)

         x = sqrt(-log(errlim)) / eta
         y = sqrt(taumax)
         ekmax = MINVAL(evr) * tpi2 * (x+y)**2 / 2d0
         call gbox(ekmax, 1d0, rb, nr1, nr2, nr3)

#ifdef DEBUG
         print *, 'Eta = ', eta
         print *, 'Previous Eta = ', etaprev
         print *, 'Recip-space sums are over ', nq1, nq2, nq3
         print *, 'Real space sums over ', nr1, nr2, nr3
#endif

         if( nq1*nq2*nq3 > 2*nr1*nr2*nr3 ) then
            if( iter == MAXIT ) then
               print *, 'Warning: possibly suboptimal value of',
     $              ' eta in Ewald sums!'
            else
               if( eta <= etaprev ) then
                  x = eta / FACT
               else
                  x = eta - (eta-etaprev) / 3d0
               end if
               etaprev = eta
               eta = x
            end if
         else if ( 2*nq1*nq2*nq3 < nr1*nr2*nr3 ) then
            if( iter == MAXIT ) then
               print *, 'Warning: possibly suboptimal value of',
     $              ' eta in Ewald sums!'
            else
               if( eta >= etaprev ) then
                  x = eta * FACT
               else
                  x = eta + (etaprev-eta) / 3d0
               end if
               etaprev = eta
               eta = x
            end if
         else
#ifdef DEBUG
            print '(1x,a,f12.4)', 'Found an optimal value of eta=', eta
#endif
            exit
         end if
      end do

#ifdef DEBUG
      print '(a,3i3)', ' Real space sums over  ', nr1, nr2, nr3
      print '(a,3i3)', ' Recip-space sums over ', nq1, nq2, nq3
#endif

      return
      end


text:
C    Copyright(C) 2014-2017 National Technology & Engineering Solutions of
C    Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains 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
C    * Neither the name of NTESS 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

C $Id: mnorm.f,v 1.1 1990/11/30 11:12:25 gdsjaar Exp $
C $Log: mnorm.f,v $
C Revision 1.1  1990/11/30 11:12:25  gdsjaar
C Initial revision
C
C
CC* FILE: [.QMESH]MNORM.FOR
CC* MODIFIED BY: TED BLACKER
CC* MODIFICATION DATE: 7/6/90
CC* MODIFICATION: COMPLETED HEADER INFORMATION
C
      SUBROUTINE MNORM (MXND, XN, YN, NXL, LLL, STDLEN)
C***********************************************************************
C
C  SUBROUTINE MNORM = FINDS THE AVERAGE LENGTH OF THOSE LINES NOT MUCH
C                     LONGER THAN THE AVERAGE
C
C***********************************************************************
C
      DIMENSION NXL (2, 3 * MXND), XN (MXND), YN (MXND)
C
      STDLEN = 0.
      NUML = 0
      S = 0.0
      DO 100 L = 1, LLL
         N1 = NXL (1, L)
         N2 = NXL (2, L)
         IF (N1 .GT. 0) THEN
            D = SQRT ((XN (N1) - XN (N2)) **2 + (YN (N1) - YN (N2)) **2)
            S = S + D
            NUML = NUML + 1
         ENDIF
  100 CONTINUE
C
      IF (NUML .LE. 0) RETURN
      TOL = 1.25 * S / FLOAT (NUML)
      NUML = 0
      S = 0.0
      DO 110 L = 1, LLL
         N1 = NXL (1, L)
         N2 = NXL (2, L)
         IF (N1 .GT. 0) THEN
            D = SQRT ((XN (N1) - XN (N2)) **2 + (YN (N1) - YN (N2)) **2)
            IF (D .LT. TOL) THEN
               S = S + D
               NUML = NUML + 1
            ENDIF
         ENDIF
  110 CONTINUE
      STDLEN = S / FLOAT (NUML)
C
      RETURN
C
      END

max_stars_repo_name: mathstuf/seacas

max_stars_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_issues_repo_name: mathstuf/seacas

max_issues_repo_licenses: BSD-3-Clause NetCDF Zlib MIT

max_forks_repo_name: mathstuf/seacas

max_forks_repo_licenses: BSD-3-Clause NetCDF Zlib MIT


text:
! Notes:
! String comparison is performed between strings of different lengths.  This behavior may not be well defined (investigate).
! Associative array implementation is used.  Therefor, old entries will be overwritten for new values.
! Currently no method is offered for deleting specific entries.

module hashtableModule
    use sha1Module
    implicit none
    private

    ! default values
    integer, parameter, private :: keysize = 32, tablesize = 32, resize_factor = 2
    real, parameter, private :: loadfactor = 0.75

    type hashnode
        character(len=keysize) :: key = ''
        class(*), allocatable :: value
        type(hashnode), pointer :: next => null()
    end type hashnode

    type hashtable
        private
        type(hashnode), dimension(:), allocatable :: table
        integer :: count = 0

    contains
        private
        procedure, public :: put => hashtable_put
        procedure :: getc => hashtable_getc
        procedure :: getr => hashtable_getr
        procedure :: geti => hashtable_geti
        generic, public :: get => getc, getr, geti
        procedure, public :: dispose => hashtable_dispose
    end type hashtable

    public :: hashtable

contains
    function hash(key)
        character(len=*) :: key
        integer :: hash
        character(len=40) :: fullhash

        fullhash = SHA1(key)
        hash = ibclr(transfer(fullhash(1:4), hash), bit_size(hash)-1)
    end function hash

    subroutine dynamicresize(this)
        class(hashtable), target :: this
        real :: thisloadfactor
        type(hashnode), dimension(:), allocatable :: temptable
        integer :: n, i, count, j
        type(hashnode), pointer :: node

        n = size(this%table)
        count = this%count
        thisloadfactor = real(count, kind(thisloadfactor)) / real(n, kind(thisloadfactor))
        if (thisloadfactor > loadfactor) then
            allocate(temptable(count))
            j = 1
            do i = 1,n
                if (allocated(this%table(i)%value)) then
                    call hashnode_write(temptable(j), this%table(i)%key, this%table(i)%value)
                    node => this%table(i)
                    do while (associated(node%next))
                        j = j + 1
                        call hashnode_write(temptable(j), node%next%key, node%next%value)
                        node => node%next
                    end do
                    j = j + 1
                end if
            end do
            call this%dispose()
            allocate(this%table(resize_factor*n))
            do i = 1,count
                call this%put(temptable(i)%key, temptable(i)%value)
            end do
            deallocate(temptable)
        end if
    end subroutine dynamicresize

    subroutine hashtable_put(this, key, value)
        class(hashtable), target :: this
        character(len=*) :: key
        class(*) :: value
        integer :: index
        type(hashnode), pointer :: node

        if (len(key) > keysize) print *, 'Warning: Key, ', key, ', is longer than storage size (', keysize, ').'

        if (.not.allocated(this%table)) then
            allocate(this%table(tablesize))
            this%count = 0
        end if

        index = mod(hash(key), size(this%table))
        if (this%table(index)%key == '' .or. this%table(index)%key == key) then
            if (allocated(this%table(index)%value)) deallocate(this%table(index)%value)
            call hashnode_write(this%table(index), key, value)
        else
            node => this%table(index)
            do while (associated(node%next) .and. node%key /= key)
                node => node%next
            end do
            if (node%key == key) then
                if (allocated(node%value)) deallocate(node%value)
                call hashnode_write(node, key, value)
            else
                allocate(node%next)
                node => node%next
                call hashnode_write(node, key, value)
                nullify(node%next)
            end if
        end if
        this%count = this%count + 1
        call dynamicresize(this)
    end subroutine hashtable_put

    subroutine hashnode_write(this, key, value)
        class(hashnode) :: this
        character(len=*) :: key
        class(*) :: value
        real :: valuer
        integer :: valuei

        this%key = key
        select type(value)
            type is (character(*))
                call hashnode_value_setc(value, this%value)
            type is (real)
                valuer = value
                allocate(this%value, source=valuer)
            type is (integer)
                valuei = value
                allocate(this%value, source=valuei)
            class default
                print *, 'Error: passed object type is not supported.'
                stop
        end select
    end subroutine hashnode_write

    subroutine hashnode_value_setc(source, dest)
        character(*) :: source
        class(*), allocatable, intent(out) :: dest

        allocate(dest, source=source)
    end subroutine hashnode_value_setc

    subroutine hashtable_getc(this, key, value)
        class(hashtable), target :: this
        character(*) :: key, value
        integer :: index
        type(hashnode), pointer :: node

        if (len(key) > keysize) print *, 'Warning: Key, ', key, ', is longer than storage size (', keysize, ').'

        index = mod(hash(key), size(this%table))
        if (.not.associated(this%table(index)%next)) then
            if (allocated(this%table(index)%value)) then
                select type(source => this%table(index)%value)
                    type is (character(*))
                        value = source
                    class default
                        print *, 'Error: returned object type is not the result type, character(*).'
                        stop
                end select
            else
                print *, 'Error: key does not exist: ', key
                stop
            end if
        else
            node => this%table(index)
            do while (node%key /= key .and. associated(node%next))
                node => node%next
            end do
            if (node%key == key .and. allocated(node%value)) then
                select type(source => this%table(index)%value)
                    type is (character(*))
                        value = source
                    class default
                        print *, 'Error: returned object type is not the result type, character(*).'
                        stop
                end select
            else
                print *, 'Error: key does not exist: ', key
                stop
            end if
        end if
    end subroutine hashtable_getc

    subroutine hashtable_getr(this, key, value)
        class(hashtable), target :: this
        character(*) :: key
        real :: value
        integer :: index
        type(hashnode), pointer :: node

        if (len(key) > keysize) print *, 'Warning: Key, ', key, ', is longer than storage size (', keysize, ').'

        index = mod(hash(key), size(this%table))
        if (.not.associated(this%table(index)%next)) then
            if (allocated(this%table(index)%value)) then
                select type(source => this%table(index)%value)
                    type is (real)
                        value = source
                    class default
                        print *, 'Error: returned object type is not the result type, real.'
                        stop
                end select
            else
                print *, 'Error: key does not exist: ', key
                stop
            end if
        else
            node => this%table(index)
            do while (node%key /= key .and. associated(node%next))
                node => node%next
            end do
            if (node%key == key .and. allocated(node%value)) then
                select type(source => this%table(index)%value)
                    type is (real)
                        value = source
                    class default
                        print *, 'Error: returned object type is not the result type, real.'
                        stop
                end select
            else
                print *, 'Error: key does not exist: ', key
                stop
            end if
        end if
    end subroutine hashtable_getr

    subroutine hashtable_geti(this, key, value)
        class(hashtable), target :: this
        character(*) :: key
        integer :: value
        integer :: index
        type(hashnode), pointer :: node

        if (len(key) > keysize) print *, 'Warning: Key, ', key, ', is longer than storage size (', keysize, ').'

        index = mod(hash(key), size(this%table))
        if (.not.associated(this%table(index)%next)) then
            if (allocated(this%table(index)%value)) then
                select type(source => this%table(index)%value)
                    type is (integer)
                        value = source
                    class default
                        print *, 'Error: returned object type is not the result type, integer.'
                        stop
                end select
            else
                print *, 'Error: key does not exist: ', key
                stop
            end if
        else
            node => this%table(index)
            do while (node%key /= key .and. associated(node%next))
                node => node%next
            end do
            if (node%key == key .and. allocated(node%value)) then
                select type(source => this%table(index)%value)
                    type is (integer)
                        value = source
                    class default
                        print *, 'Error: returned object type is not the result type, integer.'
                        stop
                end select
            else
                print *, 'Error: key does not exist: ', key
                stop
            end if
        end if
    end subroutine hashtable_geti

    subroutine hashtable_dispose(this)
        class(hashtable), target :: this
        integer :: i
        type(hashnode), pointer :: node

        if (allocated(this%table)) then
            do i = 1,size(this%table)
                node => this%table(i)
                call hashnode_clear(node)
            end do
            deallocate(this%table)
            this%count = 0
        end if
    end subroutine hashtable_dispose

    recursive subroutine hashnode_clear(this)
        type(hashnode), pointer :: this

        if (associated(this%next)) then
            call hashnode_clear(this%next)
            deallocate(this%next)
            nullify(this%next)
        end if
        if (allocated(this%value)) deallocate(this%value)
    end subroutine hashnode_clear
end module hashtableModule

max_stars_repo_name: jaketimothy/fortran-datastructures

max_issues_repo_name: jaketimothy/fortran-datastructures

max_forks_repo_name: jaketimothy/fortran-datastructures


text:
** Copyright (c) 1989, NVIDIA CORPORATION.  All rights reserved.
**
** Licensed under the Apache License, Version 2.0 (the "License");
** you may not use this file except in compliance with the License.
** You may obtain a copy of the License at
**
**     http://www.apache.org/licenses/LICENSE-2.0
**
** Unless required by applicable law or agreed to in writing, software
** distributed under the License is distributed on an "AS IS" BASIS,
** WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
** See the License for the specific language governing permissions and
** limitations under the License.

*   POINTER statements - adjustable arrays, loc intrinsic

      parameter(N = 13)
      integer result(N), expect(N)
      pointer (ptr, idum)

      data expect /
     +  2, 3, 3, 4, 			! tests 1-4, ment
     + 2, 3, 4, 3, 4, 5, 4, 5, 6	! tests 5-13, ent
     + /

      ptr = loc(result(1))
      call ment(ptr, 2)

      ptr = loc(result(5))
      call ent(ptr, 3)

      call check(result, expect, N)
      end

      subroutine ment(p, n)
      entry foo(p, n)
      dimension iarr(n, n)
      pointer (p, iarr)	! object refd by common arg declared after entry
      do i = 1, n
	  do j = 1, n
	      iarr(i, j) = i + j
	  enddo
      enddo
      return
      end
      subroutine bar(p, n)
      pointer (p, iarr)	! object refd by common arg declared before entry
      integer iarr(n, n)
      entry ent(p, n)
      do i = 1, n
	  do j = 1, n
	      iarr(i, j) = i + j
	  enddo
      enddo
      return
      end

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
      PROGRAM ALIAS
      INTEGER I,A(5)
      DO I=1,5
         A(I)=2
      ENDDO
      WRITE(*,*) 'BEFORE:',A
      CALL SUB(A,A(1))
      WRITE(*,*) 'AFTER: ' ,A
      END
      SUBROUTINE SUB(X,Y)
      INTEGER I,X(5),Y
      DO I=1,5
         X(I)=Y*X(I)
      ENDDO
      RETURN
      END


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2019 Guido Dhondt
!
!     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 2);
!
!
!     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, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine cattet(kontet,netet_,ifac,ne,ipkon,kon,ifatet,ifreetet,
     &  bc,itetfa,ifreefa,planfa,ipofa,cotet,cg,ipoeln,ieln,ifreeln,
     &  lakon,kontetor,iquad)
!
!     catalogueing the tetrahedral elements of the mesh
!
      implicit none
!
      character*8 lakon(*)
!
      integer kontet(4,*),netet_,i,ifac(4,*),ne,ipkon(*),kon(*),index,j,
     &  nodes(4),ifatet(4,*),ifreetet,itetfa(2,*),ifreefa,ipofa(*),
     &  ipoeln(*),ieln(2,*),node,ifreeln,kontetor(6,*),iquad
!
      real*8 bc(4,*),planfa(4,*),cotet(3,*),cg(3,*)
!
!     determine the first tetrahedral element or first
!     unused element
!
      do i=1,ne
         if((lakon(i)(1:4).eq.'C3D4').or.
     &      (lakon(i)(1:5).eq.'C3D10').or.
     &      (lakon(i)(1:1).eq.char(0))) exit
      enddo
      ifreetet=i
!
      do
         do j=i+1,netet_
!
!           element number supersedes largest one
!
            if(j.gt.ne) then
               kontet(4,i)=j
               i=j
               exit
            endif
!
!           tetrahedral element
!
            if((lakon(j)(1:4).eq.'C3D4').or.
     &         (lakon(j)(1:5).eq.'C3D10')) then
               kontet(4,i)=j
               i=j
               exit
            endif
!
!           unused element
!
            if(lakon(j)(1:1).eq.char(0)) then
               kontet(4,i)=j
               i=j
               exit
            endif
         enddo
         if(j.eq.netet_) exit
      enddo
      kontet(4,netet_)=0
!
!     initialization of ipofa and ifac
!
      do i=1,4*netet_
         ifac(4,i)=i+1
      enddo
      ifac(4,4*netet_)=0
!
!     initialization of ieln
!
      do i=1,4*netet_
         ieln(2,i)=i+1
      enddo
      ieln(2,4*netet_)=0
!
!     adding the tetrahedral elements one by one
!
      do i=1,ne
         if(ipkon(i).lt.0) cycle
         if((lakon(i)(1:4).ne.'C3D4').and.
     &      (lakon(i)(1:5).ne.'C3D10')) cycle
         index=ipkon(i)
         do j=1,4
            nodes(j)=kon(index+j)
         enddo
!
!        if C3D10: store the middle nodes
!        If there is at least one C3D10 element iquad is set to 1
!        which means that the refined mesh will be fully quadratic
!
         if(lakon(i)(4:4).eq.'1') then
            iquad=1
            do j=1,6
               kontetor(j,ifreetet)=kon(index+4+j)
            enddo
         endif
!
         call generatetet(kontet,ifatet,ifreetet,bc,ifac,itetfa,
     &                    ifreefa,planfa,ipofa,nodes,cotet,cg)
      enddo
!
!     generating the element per node relationship
!
      do j=1,netet_
         if(kontet(1,j).eq.0) cycle
         do i=1,4
            node=kontet(i,j)
            index=ifreeln
            ieln(1,index)=j
            ifreeln=ieln(2,index)
            if(ifreeln.eq.0) then
               write(*,*) '*ERROR in cattet: increase the'
               write(*,*) '       dimension of ieln'
               call exit(201)
            endif
            ieln(2,index)=ipoeln(node)
            ipoeln(node)=index
         enddo
      enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
      SUBROUTINE TRHT1B(IOF,DELTA)
C
C
      DOUBLE PRECISION        QBLOCK(6),          BLOCK(2), BLK(2)
C
      INTEGER                 MCB(7),   NAME(2),  IQBLK(12),IBLOCK(11)
C
      COMMON /BLANK/    BETA,     TABS,     NORAD,    RADLIN
      COMMON /TRHTX /         IK(7),    IB(7),    ICR1,     ICR2,
     1                        ICR3,     ICR4,     ICR5,     ISYM,
     2                        ICR6,     ICR7
C
      EQUIVALENCE             ( IQBLK(1),         QBLOCK(1) )
      EQUIVALENCE             ( IQBLK(2),         IBLOCK(1) )
      EQUIVALENCE             ( QBLOCK(2),        BLOCK(1) )
      EQUIVALENCE             ( QBLOCK(5),        BLK(1) )
C
C ----------------------------------------------------------------------
C
      IBLOCK(1) =2
      BLOCK(1)= 1.0D0/DELTA
      BLOCK(2)= 0.0D0
      IBLOCK(7)= 2
      BLK(1) = BETA
      BLK(2) = 0.0D0
      CALL SSG2C(IB,IK,ICR6,1,IBLOCK)
      MCB(1)=ICR6
      CALL RDTRL(MCB(1))
      IF ( MCB(4) .EQ. 6) GO TO  10
      CALL FACTRU(*40,ICR6,ICR1,ICR2,ICR3,ICR4,ICR7)
      ISYM = 0
      GO TO  20
C
C     SYMMETRIC DECOMP
C
   10 CALL FACTOR( ICR6, ICR1, ICR2, ICR3, ICR4, ICR7 )
      ISYM =1
C
C     LLL  IS ON ICR1
C
C     FORM  A  MATRIX
C
   20 BLK(1) = -(1.0D0 - BETA)
      BLK(2) = 0.0
      CALL SSG2C(IB,IK,ICR6,1,IBLOCK)
   30 RETURN
   40 CALL MESAGE(-5,ICR6,NAME)
      GO TO 30
      END


text:
C#######################################################################
C#
C#    SUBROUTINE ORTHO
C#
C#    CHECKS ORTHOGONALITY OF BASISFUNCTIONS
C#
C#    CALLED BY: STSWM
C#
      SUBROUTINE ORTHO
C
C***********************************************************************
C     DEFINE PARAMETERS
      INCLUDE 'params.i'
C***********************************************************************
C
C     GLOBAL QUANTITIES (COMMON BLOCKS)
C
C***********************************************************************
C     TRANSFORM ARRAYS
      INCLUDE 'trnsfm.i'
C
C***********************************************************************
C
C     DUMMY STORAGE (SUBROUTINE SHTRNS)
C
C     ORTHONORMALITY/ORTHOGONALITY CHECKSUMS
C
      REAL ALPON(NALP),ALPOG(NALP),TMP(0:NN)
      REAL PLOTON(0:MM,0:NN),PLOTOG(0:MM,0:NN)
C
C     TEMPORARIES
C
      INTEGER NL, JM, JN, IS, ISO
      REAL COEFF1,COEFF2, VMAX, VMIN
C
C-----------------------------------------------------------------------
C
C     INITIALIZE ORTHONORMAL ARRAY
C
      DO 100 IS = 1, NALP
         ALPON(IS) = 0.0
         ALPOG(IS) = 0.0
  100 CONTINUE
C
C     COMPUTE ORTHONORMAL SUM
C
      DO 300 NL = 1, NLAT/2
         DO 200 IS = 1, NALP
            ALPON(IS) = ALPON(IS) + ALP(IS,NL)*ALP(IS,NL)*WTS(NL)
  200    CONTINUE
  300 CONTINUE
C
      DO 390 JN = 0, NN
         IS = LDIAG(JN,2)
         DO 389 JM = 1,LDIAG(JN,1)
         IF (MOD(JN,2) .EQ. 0) THEN
            COEFF1 = 1.0
         ELSE
            COEFF1 = -1.0
         ENDIF
C
C        INITIALIZE
C
         DO 330 ISO = 0,JN-1
            TMP(ISO) = 0.0
  330    CONTINUE
C
C        INTEGRATE
C
         DO 380 NL = 1, NLAT/2
            DO 350 ISO = 0,JN-1
               IF (MOD(ISO,2) .EQ. 0) THEN
                  COEFF2 = 1.0
               ELSE
                  COEFF2 = -1.0
               ENDIF
               TMP(ISO) = TMP(ISO) +
     &         (1.0+COEFF1*COEFF2)*ALP(LDIAG(ISO,2)+JM,NL)
     &                    *ALP(IS+JM,NL)* WTS(NL)
  350       CONTINUE
  380    CONTINUE
C
C        GET MAXIMUM
C
         VMAX = 0.0
         VMIN = 0.0
         DO 385 ISO = 0, JN-1
            VMIN = MIN(TMP(ISO),VMIN)
            VMAX = MAX(TMP(ISO),VMAX)
  385    CONTINUE
         IF (ABS(VMIN) .GT. ABS(VMAX)) THEN
            ALPOG(IS+JM) = VMIN
         ELSE
            ALPOG(IS+JM) = VMAX
         ENDIF
  389    CONTINUE
  390 CONTINUE
C
C     NORMALIZE
C
      DO 400 IS = 1, NALP
         ALPON(IS) = 2.0*ALPON(IS)
  400 CONTINUE
C
C     PLOT FIELDS & OUTPUT VALUES
C
      DO 600 JN=0,NN
         DO 550 JM=0,MM
            PLOTON(JM,JN) = 0.0
            PLOTOG(JM,JN) = 0.0
  550    CONTINUE
  600 CONTINUE

      DO 500 JN = 0,NN
         IS = LDIAG(JN,2)
         WRITE (6,999) (JN+JM-1,JM-1,ALPON(IS+JM)-1.0,ALPOG(IS+JM),
     &                  JM=1,LDIAG(JN,1))
         DO 525 JM=1,LDIAG(JN,1)
            PLOTON(JM-1,JN+JM-1) = ALPON(IS+JM)-1.0
            PLOTOG(JM-1,JN+JM-1) = ALPOG(IS+JM)
  525    CONTINUE
  500 CONTINUE
C
      CALL CPCNRC(PLOTON,MM+1,MM+1,NN+1,0.0,0)
      CALL LABTOP('Orthonormality Error',0.015)
      CALL FRAME
      CALL CPCNRC(PLOTOG,MM+1,MM+1,NN+1,0.0,0)
      CALL LABTOP('Orthogonality Error',0.015)
      CALL FRAME
C
      RETURN
C
  999 FORMAT('N = ',I4,' M = ',I4,' ALPON-1.0 = ',F18.16,
     &       ' ALPOG = ',F18.16)
      END


text:
Describe Users/Sushmita here.

20100428 12:51:12 nbsp Welcome to the Wiki. Please show respect for the contributions of other editors. Users/JasonAller

20100505 07:27:59 nbsp Hey, I deleted the whole mess on Casa de Flores.  If you disagree, you can restore it by clicking Info and reverting back a couple versions... but right now theres an opportunity to start out again with a clean slate.  It seems you have some (ex?) roommate issues that spilled over onto the wiki.  Im sorry about that, but theres not much to do, as this is an open forum.  I thought this might help, considering some of the personal things that were written in public. Users/JabberWokky Evan JabberWokky Edwards


text:
	SUBROUTINE GDTXSP ( gdfile, gfunc, npts, iscale, gdtime,
     +			    gvcord, nlvl, clvl, dmin, dmax, icolor,
     +			    iline, ilwid, labflg, device, panel, gvect,
     +			    wind, havcon, havvec, iret )
C************************************************************************
C* GDTXSP								*
C*									*
C* This subroutine allows the user to exit from GDCROSS.		*
C*									*
C* GDTXSP  ( GDFILE, GFUNC, NPTS, ISCALE, GDTIME, GVCORD, NLVL,		*
C*	     CLVL, DMIN, DMAX, ICOLOR, ILINE, ILWID, LABFLG, DEVICE,	*
C*	     PANEL, GVECT, WIND, HAVCON, HAVVEC, IRET )			*
C*									*
C*									*
C* Input parameters:							*
C*	GDFILE		CHAR*		Grid file name			*
C*	GFUNC		CHAR*		Grid function			*
C*	NPTS		INTEGER		Number of points		*
C*	ISCALE		INTEGER		Scaling factor			*
C*	GDTIME		CHAR*		Grid dat/time			*
C*	GVCORD		CHAR*		Grid vertical coordinate	*
C*	NLVL		INTEGER		Number of contour levels	*
C*	CLVL  (NLVL)	REAL		Contour levels			*
C*	DMIN		REAL 		Minimum data value		*
C*	DMAX		REAL 		Maximum data value		*
C*	ICOLOR(NLVL)	INTEGER		Contour colors			*
C*	ILINE (NLVL)	INTEGER		Contour line type		*
C* 	ILWID (NLVL)	INTEGER		Contour line width		*
C*	LABFLG(NLVL)	INTEGER		Label flag			*
C*	DEVICE		CHAR*		Graphics device			*
C*	PANEL		CHAR*		Graphics panel location		*
C*	GVECT		CHAR*		Vector diagnostic function	*
C*	WIND		CHAR*		Vector plot type		*
C*	HAVCON		LOGICAL						*
C*	HAVVEC		LOGICAL						*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*					  1 = user entered EXIT		*
C*					  0 = normal return		*
C**									*
C* Log:									*
C* K. F. Brill/GSC	 6/89   Created from GDPDSP AND GDNDSP		*
C* S. Schotz/GSC	 6/90	Get respnd locally from IP_RESP		*
C************************************************************************
	CHARACTER*(*)	gdfile, gfunc, gdtime, gvcord,
     +			device, panel, gvect, wind
	INTEGER		 icolor (*), iline (*), ilwid (*), labflg (*)
	REAL		clvl (*)
	LOGICAL		respnd, havcon, havvec
C------------------------------------------------------------------------
	iret = 0
C
C*	Write information to terminal.
C
	WRITE  ( 6, 1000 )
1000	FORMAT ( ' GDTHGT PARAMETERS' / )
	CALL ST_LSTR  ( gdfile, lenf, ier )
	CALL ST_LSTR  ( gdtime, lent, ier )
	CALL ST_LSTR  ( gvcord, lenv, ier )
	CALL ST_LSTR  ( gfunc,  leng, ier )
	CALL ST_LSTR  ( device, lend, ier )
	CALL ST_LSTR  ( panel, lenp, ier )
	CALL ST_LSTR  ( gvect, lenu, ier )
	CALL ST_LSTR  ( wind, lenw, ier )
	WRITE  ( 6, 1005, IOSTAT = iostat ) gdfile ( : lenf ),
     +					    gdtime ( : lent ),
     +					    gvcord ( : lenv ),
     +					    gfunc ( : leng ),
     +					    gvect ( : lenu ),
     +					    iscale, panel, dmax, dmin
1005	FORMAT ( ' Grid file:                 ', A /
     +           ' Date/time:                 ', A /
     +           ' Vertical coordinate:       ', A /
     +           ' Grid function:             ', A /
     +           ' Vector function:           ', A /
     +           ' Scaling factor:        10**', I2 ,5x,
     +           ' Panel:  ',                    A /
     +           ' Max: ',F10.4,' -- Min: ',F10.4)
C
C*	Write out levels in groups of 7.
C
	IF ( havcon ) THEN
	    nblk = ( nlvl - 1 ) / 7 + 1
	    i1 = 1
	    i2 = 7
	    DO i = 1, nblk
		IF ( i2 .gt. nlvl ) i2 = nlvl
		WRITE ( 6, 1010 ) ( clvl   (j), j = i1, i2 )
		WRITE ( 6, 1011 ) ( icolor (j), j = i1, i2 )
		WRITE ( 6, 1012 ) ( iline  (j), j = i1, i2 )
		WRITE ( 6, 1013 ) ( ilwid  (j), j = i1, i2 )
		WRITE ( 6, 1014 ) ( labflg (j), j = i1, i2 )
1010		FORMAT ( /' LEVELS:    ',7F9.2 )
1011		FORMAT (  ' COLORS:    ',7I9 )
1012		FORMAT (  ' LINTYP:    ',7I9 )
1013		FORMAT (  ' LINWID:    ',7I9 )
1014		FORMAT (  ' LABEL:     ',7I9 )
		i1 = i1 + 7
		i2 = i2 + 7
	    END DO
	END IF
C
C*	If respond is set, wait for user to accept parameters.
C
	CALL IP_RESP ( respnd, ier )
	IF  ( respnd )  THEN
	    CALL TM_ACCP  ( ier )
	    IF  ( ier .eq. 2 )  iret = 3
	END IF
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
Describe Users/GillianGreen here.

20121226 19:00:06 nbsp Welcome to the Wiki/Business Owner Welcome to the Wiki, Gillian! Id like to suggest moving your pics to maximize the empty space in the upper right corner. If you need me to, I can help you with this. Please let me know. Users/PeteB


text:
	SUBROUTINE HSROAM  ( irmflg, ipwdth, iphght, iret )
C************************************************************************
C* HSROAM - XWP								*
C* 									*
C* This subroutine sets the roam on current driver.  			*
C* 									*
C* HSROAM  ( IRMFLG, IPWDTH, IPHGHT, IRET )				*
C*									*
C* Input parameters:                                                    *
C*      IRMFLG          INTEGER         Roam flag                       *
C*      IPWDTH          INTEGER         Pixmap(device) width            *
C*      IPHGHT          INTEGER         Pixmap(device) height           *
C*                                                                      *
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* C. Lin/EAI	 	 6/97						*
C* S. Maxwell/GSC 	 6/97	Documentation changes			*
C************************************************************************
        INCLUDE         'DEVCHR.CMN'
	INCLUDE		'ERROR.PRM'
C------------------------------------------------------------------------
        iret = NORMAL
C
C*      Check the requested device.
C*      Set roam in XW only.
C
        IF  ( ( ddev .eq. 'XW' ) .or. ( ddev .eq. 'XWP' ) )  THEN
C
            CALL XSROAM ( irmflg, ipwdth, iphght, idwdth, idhght, iret )
C
C*          Set the device size and offsets
C
            ibot   = idhght - 1
            iright = idwdth - 1
C
            isxoff = 0
            isyoff = 0
C
        END IF
C*
	RETURN
	END


text:
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!                          Futility Development Group                          !
!                             All rights reserved.                             !
!                                                                              !
! Futility is a jointly-maintained, open-source project between the University !
! of Michigan and Oak Ridge National Laboratory.  The copyright and license    !
! can be found in LICENSE.txt in the head directory of this repository.        !
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!> @brief Module provides an interface to Anderson Acceleration
!>
!> Currently supported TPLs include:
!>  - Trilinos
!>
!> @par Module Dependencies
!>
!> @par EXAMPLES
!> @code
!>
!> @endcode
!>
!> @author Ben Collins
!>   @date 07/30/2016
!>
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
MODULE trilinos_interfaces
  USE ISO_C_BINDING
#ifdef FUTILITY_HAVE_Trilinos
  USE ForTeuchos_ParameterList
#endif

  INTERFACE
!-------------------------------------------------------------------------------
! Kokkos
!-------------------------------------------------------------------------------
#ifdef FUTILITY_HAVE_Kokkos
    SUBROUTINE InitializeKokkos(num_threads,device_id) bind(C,NAME="InitializeKokkos")
      USE ISO_C_BINDING
      INTEGER(C_INT),INTENT(IN),VALUE :: num_threads
      INTEGER(C_INT),INTENT(IN),VALUE :: device_id
    ENDSUBROUTINE

    SUBROUTINE FinalizeKokkos() bind(C,NAME="FinalizeKokkos")
      USE ISO_C_BINDING
    ENDSUBROUTINE
#endif

!-------------------------------------------------------------------------------
! Vector Interfaces
!-------------------------------------------------------------------------------
    SUBROUTINE ForPETRA_VecInit(id,n,nlocal,comm) bind(C,NAME="ForPETRA_VecInit")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(INOUT)    :: id
      INTEGER(C_INT),INTENT(IN),VALUE :: n
      INTEGER(C_INT),INTENT(IN),VALUE :: nlocal
      INTEGER(C_INT),INTENT(IN),VALUE :: comm
    ENDSUBROUTINE

    SUBROUTINE ForPETRA_VecSetImportMap(id,n,gids) bind(C,NAME="ForPETRA_VecSetImportMap")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(IN),VALUE        :: id
      INTEGER(C_INT),INTENT(IN),VALUE        :: n
      INTEGER(C_INT),DIMENSION(*),INTENT(IN) :: gids
    ENDSUBROUTINE

    SUBROUTINE ForPETRA_VecDestroy(id) bind(C,NAME="ForPETRA_VecDestroy")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(IN),VALUE :: id
    ENDSUBROUTINE

    SUBROUTINE ForPETRA_VecSet(id,i,val) bind(C,NAME="ForPETRA_VecSet")
      IMPORT :: C_INT,C_DOUBLE
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      INTEGER(C_INT),INTENT(IN),VALUE  :: i
      REAL(C_DOUBLE),INTENT(IN),VALUE  :: val
    ENDSUBROUTINE

    SUBROUTINE ForPETRA_VecSetAll(id,val) bind(C,NAME="ForPETRA_VecSetAll")
      IMPORT :: C_INT,C_DOUBLE
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      REAL(C_DOUBLE),INTENT(IN),VALUE  :: val
    ENDSUBROUTINE

    SUBROUTINE ForPETRA_VecTransfer(id) bind(C,NAME="ForPETRA_VecTransfer")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(IN),VALUE :: id
    ENDSUBROUTINE

    SUBROUTINE ForPETRA_VecGet(id,i,val) bind(C,NAME="ForPETRA_VecGet")
      IMPORT :: C_INT,C_DOUBLE
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      INTEGER(C_INT),INTENT(IN),VALUE  :: i
      REAL(C_DOUBLE),INTENT(OUT)       :: val
    ENDSUBROUTINE

    SUBROUTINE ForPETRA_VecCopy(id,idfrom) bind(C,NAME="ForPETRA_VecCopy")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      INTEGER(C_INT),INTENT(IN),VALUE  :: idfrom
    ENDSUBROUTINE

    SUBROUTINE ForPETRA_VecAXPY(id,idx,a,b) bind(C,NAME="ForPETRA_VecAXPY")
      IMPORT :: C_INT,C_DOUBLE
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      INTEGER(C_INT),INTENT(IN),VALUE  :: idx
      REAL(C_DOUBLE),INTENT(IN),VALUE  :: a
      REAL(C_DOUBLE),INTENT(IN),VALUE  :: b
    ENDSUBROUTINE

    SUBROUTINE ForPETRA_VecSum(id,val) bind(C,NAME="ForPETRA_VecSum")
      IMPORT :: C_INT,C_DOUBLE
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      REAL(C_DOUBLE),INTENT(OUT)       :: val
    ENDSUBROUTINE

    SUBROUTINE ForPETRA_VecNorm2(id,val) bind(C,NAME="ForPETRA_VecNorm2")
      IMPORT :: C_INT,C_DOUBLE
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      REAL(C_DOUBLE),INTENT(OUT)       :: val
    ENDSUBROUTINE

    SUBROUTINE ForPETRA_VecMax(id,val) bind(C,NAME="ForPETRA_VecMax")
      IMPORT :: C_INT,C_DOUBLE
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      REAL(C_DOUBLE),INTENT(OUT)       :: val
    ENDSUBROUTINE

    SUBROUTINE ForPETRA_VecScale(id,val) bind(C,NAME="ForPETRA_VecScale")
      IMPORT :: C_INT,C_DOUBLE
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      REAL(C_DOUBLE),INTENT(IN),VALUE  :: val
    ENDSUBROUTINE

    ! When calling from FORTRAN, need to append C_NULL_CHAR to name:
    !    "myfile.out" // c_null_char
    SUBROUTINE ForPETRA_VecEdit(id,name) bind(C,NAME="ForPETRA_VecEdit")
      IMPORT :: C_INT, C_CHAR
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      CHARACTER(KIND=C_CHAR),INTENT(IN),DIMENSION(*) :: name
    ENDSUBROUTINE

!-------------------------------------------------------------------------------
! Matrix Interfaces
!-------------------------------------------------------------------------------
    SUBROUTINE ForPETRA_MatInit(id,n,nlocal,rnnz,comm) bind(C,NAME="ForPETRA_MatInit")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(INOUT)    :: id
      INTEGER(C_INT),INTENT(IN),VALUE :: n
      INTEGER(C_INT),INTENT(IN),VALUE :: nlocal
      INTEGER(C_INT),INTENT(IN),VALUE :: rnnz
      INTEGER(C_INT),INTENT(IN),VALUE :: comm
    ENDSUBROUTINE

    SUBROUTINE ForPETRA_MatDestroy(id) bind(C,NAME="ForPETRA_MatDestroy")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(IN),VALUE        :: id
    ENDSUBROUTINE

    SUBROUTINE ForPETRA_MatReset(id) bind(C,NAME="ForPETRA_MatReset")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(IN),VALUE        :: id
    ENDSUBROUTINE

    SUBROUTINE ForPETRA_MatSet(id,i,nnz,j,val) bind(C,NAME="ForPETRA_MatSet")
      IMPORT :: C_INT,C_DOUBLE
      INTEGER(C_INT),INTENT(IN),VALUE        :: id
      INTEGER(C_INT),INTENT(IN),VALUE        :: i
      INTEGER(C_INT),INTENT(IN),VALUE        :: nnz
      INTEGER(C_INT),DIMENSION(*),INTENT(IN) :: j
      REAL(C_DOUBLE),DIMENSION(*),INTENT(IN) :: val
    ENDSUBROUTINE

    SUBROUTINE ForPETRA_MatGet(id,i,j,val) bind(C,NAME="ForPETRA_MatGet")
      IMPORT :: C_INT,C_DOUBLE
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      INTEGER(C_INT),INTENT(IN),VALUE  :: i
      INTEGER(C_INT),INTENT(IN),VALUE  :: j
      REAL(C_DOUBLE),INTENT(OUT)       :: val
    ENDSUBROUTINE

    SUBROUTINE ForPETRA_MatAssemble(id) bind(C,NAME="ForPETRA_MatAssemble")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
    ENDSUBROUTINE

    SUBROUTINE ForPETRA_MatMult(idA,trans,idX,idY) bind(C,NAME="ForPETRA_MatMult")
      IMPORT :: C_INT,C_BOOL
      INTEGER(C_INT),INTENT(IN),VALUE        :: idA
      LOGICAL(C_BOOL),INTENT(IN),VALUE       :: trans
      INTEGER(C_INT),INTENT(IN),VALUE        :: idX
      INTEGER(C_INT),INTENT(IN),VALUE        :: idY
    ENDSUBROUTINE

    ! When calling from FORTRAN, need to append C_NULL_CHAR to name:
    !    "myfile.out" // c_null_char
    SUBROUTINE ForPETRA_MatEdit(id,name) bind(C,NAME="ForPETRA_MatEdit")
      IMPORT :: C_INT,C_CHAR
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      CHARACTER(KIND=C_CHAR),INTENT(IN),DIMENSION(*) :: name
    ENDSUBROUTINE

!-------------------------------------------------------------------------------
! Anasazi Interfaces
!-------------------------------------------------------------------------------
    SUBROUTINE Anasazi_Init(id) bind(C,NAME="Anasazi_Init")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(INOUT)    :: id
    ENDSUBROUTINE

#ifdef FUTILITY_HAVE_Trilinos
    SUBROUTINE Anasazi_Init_Params(id, plID) bind(C,NAME="Anasazi_Init_Params")
      IMPORT :: C_INT
      IMPORT :: ForTeuchos_ParameterList_ID
      INTEGER(C_INT),INTENT(INOUT)    :: id
      TYPE(ForTeuchos_ParameterList_ID), INTENT(IN) :: plID
    ENDSUBROUTINE
#endif

    SUBROUTINE Anasazi_Destroy(id) bind(C,NAME="Anasazi_Destroy")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
    ENDSUBROUTINE

    SUBROUTINE Anasazi_SetMat(id,idLHS,idRHS) bind(C,NAME="Anasazi_SetMat")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      INTEGER(C_INT),INTENT(IN),VALUE  :: idLHS
      INTEGER(C_INT),INTENT(IN),VALUE  :: idRHS
    ENDSUBROUTINE

    SUBROUTINE Anasazi_SetPC(id,idPC) bind(C,NAME="Anasazi_SetPC")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      INTEGER(C_INT),INTENT(IN),VALUE  :: idPC
    ENDSUBROUTINE

    SUBROUTINE Anasazi_SetX(id,idX) bind(C,NAME="Anasazi_SetX")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      INTEGER(C_INT),INTENT(IN),VALUE  :: idX
    ENDSUBROUTINE

    SUBROUTINE Anasazi_SetConvCrit(id,tol,maxit) bind(C,NAME="Anasazi_SetConvCrit")
      IMPORT :: C_INT,C_DOUBLE
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      REAL(C_DOUBLE),INTENT(IN),VALUE  :: tol
      INTEGER(C_INT),INTENT(IN),VALUE  :: maxit
    ENDSUBROUTINE

    SUBROUTINE Anasazi_Solve(id) bind(C,NAME="Anasazi_Solve")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
    ENDSUBROUTINE

    SUBROUTINE Anasazi_GetEigenvalue(id,k) bind(C,NAME="Anasazi_GetEigenvalue")
      IMPORT :: C_INT,C_DOUBLE
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      REAL(C_DOUBLE),INTENT(OUT)       :: k
    ENDSUBROUTINE

    SUBROUTINE Anasazi_GetResid(id,resid) bind(C,NAME="Anasazi_GetResid")
      IMPORT :: C_INT,C_DOUBLE
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      REAL(C_DOUBLE),INTENT(OUT)       :: resid
    ENDSUBROUTINE

    SUBROUTINE Anasazi_GetIterationCount(id,niter) bind(C,NAME="Anasazi_GetIterationCount")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      INTEGER(C_INT),INTENT(OUT)       :: niter
    ENDSUBROUTINE

!-------------------------------------------------------------------------------
! Belos Interfaces
!-------------------------------------------------------------------------------
    SUBROUTINE Belos_Init(id) bind(C,NAME="Belos_Init")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(INOUT)    :: id
    ENDSUBROUTINE

#ifdef FUTILITY_HAVE_Trilinos
    SUBROUTINE Belos_Init_Params(id, plID) bind(C,NAME="Belos_Init_Params")
      IMPORT :: C_INT
      IMPORT :: ForTeuchos_ParameterList_ID
      INTEGER(C_INT),INTENT(INOUT)    :: id
      TYPE(ForTeuchos_ParameterList_ID), INTENT(IN) :: plID
    ENDSUBROUTINE
#endif

    SUBROUTINE Belos_Destroy(id) bind(C,NAME="Belos_Destroy")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
    ENDSUBROUTINE

    SUBROUTINE Belos_SetMat(id,idA) bind(C,NAME="Belos_SetMat")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      INTEGER(C_INT),INTENT(IN),VALUE  :: idA
    ENDSUBROUTINE

    SUBROUTINE Belos_SetPC(id,idPC) bind(C,NAME="Belos_SetPC")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      INTEGER(C_INT),INTENT(IN),VALUE  :: idPC
    ENDSUBROUTINE

    SUBROUTINE Belos_SetX(id,idX) bind(C,NAME="Belos_SetX")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      INTEGER(C_INT),INTENT(IN),VALUE  :: idX
    ENDSUBROUTINE

    SUBROUTINE Belos_Setb(id,idb) bind(C,NAME="Belos_Setb")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      INTEGER(C_INT),INTENT(IN),VALUE  :: idb
    ENDSUBROUTINE

    SUBROUTINE Belos_SetConvCrit(id,tol,maxit) bind(C,NAME="Belos_SetConvCrit")
      IMPORT :: C_INT,C_DOUBLE
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      REAL(C_DOUBLE),INTENT(IN),VALUE  :: tol
      INTEGER(C_INT),INTENT(IN),VALUE  :: maxit
    ENDSUBROUTINE

    SUBROUTINE Belos_Solve(id) bind(C,NAME="Belos_Solve")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
    ENDSUBROUTINE

    SUBROUTINE Belos_GetResid(id,resid) bind(C,NAME="Belos_GetResid")
      IMPORT :: C_INT,C_DOUBLE
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      REAL(C_DOUBLE),INTENT(OUT)       :: resid
    ENDSUBROUTINE

    SUBROUTINE Belos_GetIterationCount(id,niter) bind(C,NAME="Belos_GetIterationCount")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      INTEGER(C_INT),INTENT(OUT)       :: niter
    ENDSUBROUTINE

!-------------------------------------------------------------------------------
! Preconditioner Interfaces
!-------------------------------------------------------------------------------
    SUBROUTINE Preconditioner_Init(id,opt) bind(C,NAME="Preconditioner_Init")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(INOUT)    :: id
      INTEGER(C_INT),INTENT(IN),VALUE :: opt
    ENDSUBROUTINE

#ifdef FUTILITY_HAVE_Trilinos
    SUBROUTINE Preconditioner_InitParams(id,plist) bind(C,NAME="Preconditioner_InitParams")
      IMPORT :: C_INT
      IMPORT :: ForTeuchos_ParameterList_ID
      INTEGER(C_INT),INTENT(INOUT)    :: id
      TYPE(ForTeuchos_ParameterList_ID), INTENT(IN) :: plist
    ENDSUBROUTINE
#endif

    SUBROUTINE Preconditioner_Destroy(id) bind(C,NAME="Preconditioner_Destroy")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(IN),VALUE :: id
    ENDSUBROUTINE

    SUBROUTINE Preconditioner_Setup(id,idM) bind(C,NAME="Preconditioner_Setup")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      INTEGER(C_INT),INTENT(IN),VALUE  :: idM
    ENDSUBROUTINE

    SUBROUTINE Preconditioner_Reset(id,idM) bind(C,NAME="Preconditioner_Reset")
      IMPORT :: C_INT
      INTEGER(C_INT),INTENT(IN),VALUE  :: id
      INTEGER(C_INT),INTENT(IN),VALUE  :: idM
    ENDSUBROUTINE

!-------------------------------------------------------------------------------
! Anderson Interfaces
!-------------------------------------------------------------------------------
!    SUBROUTINE Anderson_Init(id,depth,beta,start,idv) bind(C,NAME="Anderson_Init")
!      IMPORT :: C_INT,C_DOUBLE
!      INTEGER(C_INT),INTENT(INOUT)    :: id
!      INTEGER(C_INT),INTENT(IN),VALUE :: depth
!      REAL(C_DOUBLE),INTENT(IN),VALUE :: beta
!      INTEGER(C_INT),INTENT(IN),VALUE :: start
!      INTEGER(C_INT),INTENT(IN),VALUE :: idv
!    ENDSUBROUTINE
!
!    SUBROUTINE Anderson_Destroy(id) bind(C,NAME="Anderson_Destroy")
!      IMPORT :: C_INT
!      INTEGER(C_INT),INTENT(IN),VALUE :: id
!    ENDSUBROUTINE
!
!    SUBROUTINE Anderson_Update(id) bind(C,NAME="Anderson_Update")
!      IMPORT :: C_INT
!      INTEGER(C_INT),INTENT(IN),VALUE  :: id
!    ENDSUBROUTINE
!
!    SUBROUTINE Anderson_Reset(id) bind(C,NAME="Anderson_Reset")
!      IMPORT :: C_INT
!      INTEGER(C_INT),INTENT(IN),VALUE  :: id
!    ENDSUBROUTINE

!-------------------------------------------------------------------------------
! Rythmos Interfaces
!-------------------------------------------------------------------------------
!    SUBROUTINE TS_Init(id,funptr,n,tol) bind(C,NAME="TS_Init")
!      IMPORT :: C_INT,C_DOUBLE,C_FUNPTR
!      INTEGER(C_INT),INTENT(INOUT)    :: id
!      TYPE(C_FUNPTR),INTENT(IN),VALUE :: funptr
!      INTEGER(C_INT),INTENT(IN),VALUE :: n
!      REAL(C_DOUBLE),INTENT(IN),VALUE :: tol
!    ENDSUBROUTINE
!
!#ifdef FUTILITY_HAVE_Trilinos
!    SUBROUTINE TS_InitParams(id,funptr,n,tol,plist) bind(C,NAME="TS_InitParams")
!      IMPORT :: C_INT,C_DOUBLE,C_FUNPTR
!      IMPORT :: ForTeuchos_ParameterList_ID
!      INTEGER(C_INT),INTENT(INOUT)    :: id
!      TYPE(C_FUNPTR),INTENT(IN),VALUE :: funptr
!      INTEGER(C_INT),INTENT(IN),VALUE :: n
!      REAL(C_DOUBLE),INTENT(IN),VALUE :: tol
!      TYPE(ForTeuchos_ParameterList_ID), INTENT(IN) :: plist
!    ENDSUBROUTINE
!#endif
!
!    SUBROUTINE TS_Destroy(id) bind(C,NAME="TS_Destroy")
!      IMPORT :: C_INT
!      INTEGER(C_INT),INTENT(IN),VALUE :: id
!    ENDSUBROUTINE
!
!    SUBROUTINE TS_Step(id,tstart,tend,xstart,xend) bind(C,NAME="TS_Step")
!      IMPORT :: C_INT,C_DOUBLE
!      INTEGER(C_INT),INTENT(IN),VALUE  :: id
!      REAL(C_DOUBLE),INTENT(IN) :: tstart
!      REAL(C_DOUBLE),INTENT(IN) :: tend
!      REAL(C_DOUBLE),DIMENSION(*),INTENT(IN) :: xstart
!      REAL(C_DOUBLE),DIMENSION(*),INTENT(OUT) :: xend
!    ENDSUBROUTINE

  ENDINTERFACE
ENDMODULE trilinos_interfaces

max_stars_repo_name: picmc/Futility

max_issues_repo_name: picmc/Futility

max_forks_repo_name: wgurecky/Futility


text:
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!                                                                      %
! Copyright (C) 2015 Adrian Martinez Vargas                            %
!                                                                      %
! This software may be modified and distributed under the terms        %
! of the MIT license.  See the LICENSE.txt file for details.           %
!                                                                      %
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!                                                                      %
! Copyright (C) 1996, The Board of Trustees of the Leland Stanford     %
! Junior University.  All rights reserved.                             %
!                                                                      %
! The programs in GSLIB are distributed in the hope that they will be  %
! useful, but WITHOUT ANY WARRANTY.  No author or distributor accepts  %
! responsibility to anyone for the consequences of using them or for   %
! whether they serve any particular purpose or work at all, unless he  %
! says so in writing.  Everyone is granted permission to copy, modify  %
! and redistribute the programs in GSLIB, but only under the condition %
! that this notice and the above copyright notice remain intact.       %
!                                                                      %
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!-----------------------------------------------------------------------
!
! The functions and subroutines below were modified from the
! version 2.0 of the gslib code written in fortran 77
!
! The objective is to add functionality to GSLIB and to link
! this code with python using f2py. It uses f90 code convention.
! The arrays are dynamic and externally declared (from python)
!
! the code was converted from Fortran 77 to Fortran 90 using F2F.pl
!
! for more information please refer to:
! - gslib77 source code: http://www.statios.com/software/gslib77_ls.tar.gz
! - GSLIB: Geostatistical Software Library and User's Guide. Second edition
!   by Clayton V. Deutsch, Andre G. Journel, 1997.
! - F2PY Users Guide: http://docs.scipy.org/doc/numpy-dev/f2py/
! - F2F https://bitbucket.org/lemonlab/f2f
!-----------------------------------------------------------------------

!*********************************************************************************
!     Subroutines in GSLIB (auxiliary functions)
!*********************************************************************************
! include 'setrot.f90'
! include 'cova3.f90'
! include 'sqdist.f90'
! include 'getindx.f90'
include 'locate.f90'
! this includes sortem and dsortem
include 'sortem.f90'
include 'acorni.f90'
! include 'gauinv.f90'
include 'powint.f90'
! include 'gcum.f90'
include 'getz.f90'

subroutine trans(ivtype,nr,vr,wt,no,vo,wo,nx,ny,nz,wx,wy,wz, &
                 nxyza,zmin,zmax,ltpar,utpar,ltail,utail,ldata,kv,ef,rseed, &
                 gmedian,rvr,rcdf,ncut,zval,error)

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

    !                      Univariate Transformation
    !                      *************************

    ! Reads in a reference distribution and a number of other distributions
    ! and then transforms the values in each of the second distributions
    ! such that their histograms match that of the reference distribution.

    ! INPUT/OUTPUT Parameters:

    !   ivtype      variable type (1=continuous, 0=categorical)
    !   nx, ny, nz  size of categorical variable realizations to transform
    !   wx, wy, wz  window size for breaking ties
    !   nxyz        size to of continuous variable data set to transform
    !   zmin,zmax   minimum and maximum data values
    !   ltail,ltpar lower tail: option, parameter
    !   utail,utpar upper tail: option, parameter
    !   ldata       honor local data (1=yes, 0=no)
    !   wtfac       control parameter

    ! This version use only one ralization

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

    implicit none

    ! input
    integer, intent(in) :: ivtype              ! type of distribution 1=continuous, 0=categorical
    integer, intent(in) :: nr
    real, intent(in), dimension(nr) :: vr, wt  ! reference distribution, variable and weight
    integer, intent(in) :: no
    real, intent(in), dimension(no) :: vo, wo  ! original distributions variable and weight
    integer, intent(in) :: nx, ny, nz          ! categorical: nx, ny, nz: size of 3-D model
    integer, intent(in) :: wx, wy, wz             !     window size for tie-breaking
    integer, intent(in) :: nxyza                ! continuous: number to transform per "set"
    real, intent(in) :: zmin,zmax              !     minimum and maximum values
    real, intent(in) :: ltpar,utpar            !     lower, upper tail parameter
    integer, intent(in) :: ltail,utail         !     lower, upper tail function type
    integer, intent(in) :: ldata               ! honor local data? (1=yes, 0=no)
    real, intent(in), dimension(nr) :: kv      ! estimation variance
    real, intent(in) ::    ef               ! control parameter or scaling factor ( 0.33 < w < 3.0 )
    integer, intent(in) ::     rseed           ! random number seed (conditioning cat.)

    ! output
    integer, intent(out) :: error
    real, intent(out) :: gmedian                  ! median (continous)
    real, intent(out), dimension(nr) :: rvr,rcdf  ! category and cumulative proportion (categorical)
    integer, intent(out) ::                ncut      ! num of categories (use rvr[1:ncut] or rcdf[1:ncut] )
    real, intent(out), dimension(no) :: zval      ! the transformed variable

    ! internal
    real ::  dcdf(no),dvr(no), &
             indx(no),fuzzcat(no),catcdf(nr)
    integer ::   category(nr), KORDEI, MAXOP1, i, j, &
                 icut,iix,iiy,iiz,isim,iwt,iwtd,ix,iy,iz, &
                 jx,jy,jz, nloc, num, nxy, nxyz
    real :: EPSLON, wtfac, wtt, tcdf, evmax, cp, oldcp ,rn,tmax,powint,getz
    real*8 ::  acorni, ixv(12)
    real, dimension(1) :: c,d,e,f,g,h

    wtfac = ef
    EPSLON=1.0e-12
    tmax = 1.0e21
    nxyz = nxyza

    ! ACORN parameters:
    KORDEI=12 ! this is hard coded... ixv (KORDEI=12)
    MAXOP1 = KORDEI+1

    if(ldata == 1) then
        ! Scale from 0 to 1  --> 0.33 to 3.0
        wtfac = 0.33 + wtfac*(3.0-0.33)

        do i=1,MAXOP1
            ixv(i) = 0.0
        end do
        ixv(1) = rseed
        do i=1,10000
            rn = real(acorni(ixv, KORDEI))
        end do
    end if

    ! Check for error situation:

    if(ltail /= 1 .AND. ltail /= 2) then
        error =10 ! 'ERROR invalid lower tail option, only allow 1 or 2 - see manual '
        stop
    endif
    if(utail /= 1 .AND. utail /= 2 .AND. utail /= 4) then
        error =11 ! 'ERROR invalid upper tail option, only allow 1,2 or 4 - see manual '
        stop
    endif
    if(utail == 4 .AND. utpar < 1.0) then
        error =12 ! 'ERROR invalid power for hyperbolic tail must be greater than 1.0!'
        stop
    endif

    ! Read as much data for reference distribution as possible:

    ncut = 0
    tcdf = 0
    do j=1, nr

        wtt = wt(j)

        if(ivtype == 1) then
            ncut = ncut + 1
            rvr(ncut)  = vr(j)
            rcdf(ncut) = wtt
        else
            do i=1,ncut
                if(int(vr(j)+0.5) == int(rvr(i))) then
                    icut = i
                    go to 4
                end if
            end do
            ncut = ncut + 1
            rvr(ncut)  = vr(j)
            rcdf(ncut) = 0.0
            icut       = ncut
            4 continue
            rcdf(icut) = rcdf(icut) + wtt
        end if
        tcdf = tcdf + wtt

    end do

    ! Sort the Reference Distribution and Check for error situation:

    call sortem(1,ncut,rvr,1,rcdf,c,d,e,f,g,h)
    if(ncut <= 1 .OR. tcdf <= EPSLON) then
        error =100 ! 'ERROR: too few data or too low weight'
        return
    endif
    if(ivtype == 1 .AND. utail == 4 .AND. rvr(ncut) <= 0.0) then
        error = 200 !'ERROR can not use hyperbolic tail with negative values! - see manual '
        return
    endif

    ! Turn the (possibly weighted) distribution into a cdf that is useful:

    tcdf  = 1.0 / tcdf
    if(ivtype == 1) then
        oldcp = 0.0
        cp    = 0.0
        do i=1,ncut
            cp     = cp + rcdf(i) * tcdf
            rcdf(i) =(cp + oldcp) * 0.5
            oldcp  = cp
        end do
    else
        do i=1,ncut
            rcdf(i) = rcdf(i) * tcdf
        end do
    end if

    ! Write Some of the Statistics to the screen:

    if(ivtype == 1) then
        call locate(rcdf,ncut,1,ncut,0.5,j)
        gmedian = powint(rcdf(j),rcdf(j+1),rvr(j),rvr(j+1),0.5,1.0)
    else
        do i=1,ncut
            catcdf(i) = 0.0
            if(i > 1) rcdf(i) = rcdf(i) + rcdf(i-1)
        end do
    end if

    ! MAIN Loop over all the increments to transform:

    if(ivtype == 0) then
        nxyz = nx*ny*nz
        nxy  = nx*ny
    end if

    if(nxyz/=nxyza) then
        error =-1 ! 'Warning:  nx*ny*nz/=nxyza'
    endif

    if(nxyz/=nr) then
        error =-2 ! 'Warning:  nx*ny*nz/=nr'
    endif

    ! Read in the data values:

    tcdf = 0.0
    num  = 0
    do i=1,no
        num = num + 1
        if (no<num) then
            error = 1000  ! bad parameters, exceding array length no<num
            return
        end if
        dvr(num)  = vo(num)
        indx(num) = real(num)
        wtt = wo(num)
        dcdf(num) = wtt
        tcdf      = tcdf + wtt

        ! Keep track of the proportions if working with a categorical variable:

        if(ivtype == 0) then
            do j=1,ncut
                if(int(dvr(num)+0.5) == int(rvr(j))) then
                    icut = j
                    go to 6
                end if
            end do
            error=2000 ! 'Found a code not found in reference'   write(*,*) dvr(num)
            return
            6 continue
            catcdf(icut)  = catcdf(icut) + wtt
        end if
    end do

    if(tcdf <= EPSLON) then
        error=3000 ! 'ERROR: no data'
        return
    endif
    if(ivtype == 0 .AND. num /= nxyz) then
        error=400 ! 'ERROR: you must have nxyz for a categorical transformation'
        return
    endif

    ! For categorical transformation we need a fuzzy category which is the
    ! local average category - this is for breaking ties in the
    ! transformation procedure.

    if(ivtype == 0) then
        do i=1,nxyz
            fuzzcat(i) = 0.0
            nloc       = 0
            iz = int((i-1)/nxy) + 1
            iy = int((i-(iz-1)*nxy-1)/nx) + 1
            ix = i- (iz-1)*nxy - (iy-1)*nx
            do iix=-wx,wx
                do iiy=-wy,wy
                    do iiz=-wz,wz
                        jx = ix + iix
                        jy = iy + iiy
                        jz = iz + iiz
                        if(jx >= 1 .AND. jx <= nx .AND. &
                        jy >= 1 .AND. jy <= ny .AND. &
                        jz >= 1 .AND. jz <= nz) then
                            j = jx + (jy-1)*nx + (jz-1)*nxy
                            nloc = nloc + 1
                            fuzzcat(i) = fuzzcat(i) + dvr(j)
                        end if
                    end do
                end do
            end do
            if(nloc > 0) then
                fuzzcat(i) = fuzzcat(i) / real(nloc)
            else
                fuzzcat(i) = 2*tmax
            end if
        end do
    end if

    ! Turn the (possibly weighted) data distribution into a useful cdf:

    if(ivtype == 1) then
        call sortem(1,num,dvr,2,dcdf,indx,d,e,f,g,h)
    else
        call sortem(1,num,fuzzcat,3,dcdf,dvr,indx,e,f,g,h)
    end if
    oldcp = 0.0
    cp    = 0.0
    tcdf  = 1.0 / tcdf
    do i=1,num
        cp     = cp + dcdf(i) * tcdf
        dcdf(i) =(cp + oldcp) * 0.5
        oldcp  = cp
    end do

    ! Now, get the right order back:

    if(ivtype == 1) then
        call sortem(1,num,indx,2,dcdf,dvr,d,e,f,g,h)
    else
        call sortem(1,num,indx,3,dcdf,dvr,fuzzcat,e,f,g,h)
    end if

    ! Read in the kriging variance to array "indx" if we have to honor
    ! local data:

    if(ldata == 1) then
        evmax = -1.0e21
        do i=1,num
            indx(i) = kv(i)
            if(indx(i) >= 0.0) then
                indx(i) = sqrt(max(indx(i),0.0))
                if(indx(i) > evmax) evmax = indx(i)
            end if
        end do
    end if

    ! Go through all the data back transforming them to the reference CDF:

    do i=1,num

        if(ivtype == 1) then
            zval(i) = getz(dcdf(i),ncut,rvr,rcdf,zmin, &
            zmax,ltail,ltpar,utail,utpar)
        else
            zval(i) = rvr(1)
            do j=2,ncut
                if(dcdf(i) > rcdf(j-1)) &
                zval(i) = rvr(j)
            end do
        end if

        ! Now, do we have to honor local data?

        if(ldata == 1) then
            if(indx(i) >= 0) then
                wtt = (indx(i)/evmax)**wtfac
                if(ivtype == 1) then
                    zval(i) = dvr(i)+wtt*(zval(i)-dvr(i))
                else
                    rn = real(acorni(ixv, KORDEI))
                    if(rn > wtt) zval(i) = dvr(i)
                end if
            end if
        end if

    end do

end subroutine trans

max_stars_repo_name: cmrajan/pygslib

max_issues_repo_name: kaufmanno/pygslib

max_forks_repo_name: kaufmanno/pygslib