text:
C------------------------------------------------------------------------------
        subroutine ftxmsg(action,text)

C       get, put, or clear the error message stack

        integer action
        character*(*) text

        integer nbuff,i
        parameter (nbuff=50)
        character*80 txbuff(nbuff)
        save txbuff
        data txbuff/nbuff * ' '/

        if (action .eq. -1)then

C           get error message from top of stack and shift the stack up one
            text=txbuff(1)
            do 10 i=1,nbuff-1
                txbuff(i) = txbuff(i+1)
 10         continue
            txbuff(nbuff)=' '

        else if (action .eq. 1)then

C           put error message onto stack.
            do 20 i=1,nbuff
                if (txbuff(i) .eq. ' ')then
                   txbuff(i)=text
                   return
                end if
20          continue
C           stack is full so discard oldest message
            do 25 i=1,nbuff-1
                txbuff(i) = txbuff(i+1)
25          continue
            txbuff(nbuff)=text

        else if (action .eq. 0)then

C           clear the error message stack
            do 30 i=1,nbuff
                txbuff(i) = ' '
30          continue

        end if
        end

ext: f

max_stars_repo_name: iraf-community/stsdas

max_issues_repo_name: kirxkirx/iraf

max_forks_repo_name: kirxkirx/iraf


!--------------------------------------------
! Unified Pseudopotential Format (UPF)
! This potential is adopted in QUANTUM ESPRESSO
! The unit of energy is in Rydberg, and converted to Hartree in this routine
!--------------------------------------------
module read_paw_upf_module

  use structures,only : s_pp_info
  use paw_variables,only : paw
  use radial_1d
  use builtin_pz, only: PZxc

  implicit none

  private
  public :: read_paw_upf
  public :: get_mr_norb_upf

contains

  subroutine read_paw_upf( g, pp, ik )
    implicit none
    integer,intent(in) :: g, ik
    type(s_pp_info),intent(inout) :: pp
    character(200) :: cbuf
    integer :: nrr, norb, nelem, i, j, l, n, ll,l0
    integer :: nwfc
    real(8),allocatable :: work(:),qij(:)
    real(8),allocatable :: vh(:),vxc(:)
    real(8) :: tmp,exc,trho,pi4,tmp2,tmp1

    nrr=maxval( pp%mr )
    norb=maxval( pp%num_orb )
    nelem=size( pp%mr )

    write(*,*) "nrr,norb,nelem",nrr,norb,nelem
    allocate( qij(nrr) ); qij=0.0d0

    allocate( paw%rad(nrr,nelem)       ); paw%rad=0.0d0
    allocate( paw%rab(nrr,nelem)       ); paw%rab=0.0d0
    allocate( paw%aewf(nrr,norb,nelem) ); paw%aewf=0.0d0
    allocate( paw%pswf(nrr,norb,nelem) ); paw%pswf=0.0d0
    allocate( paw%beta(nrr,norb,nelem) ); paw%beta=0.0d0
    allocate( paw%occ(norb,nelem)      ); paw%occ=0.0d0
    allocate( paw%rhops(nrr,nelem)     ); paw%rhops=0.0d0
    allocate( paw%rhoae(nrr,nelem)     ); paw%rhoae=0.0d0
    allocate( paw%rhocps(nrr,nelem)    ); paw%rhocps=0.0d0
    allocate( paw%rhocae(nrr,nelem)    ); paw%rhocae=0.0d0
    allocate( paw%vlocps(nrr,nelem)    ); paw%vlocps=0.0d0
    allocate( paw%vlocae(nrr,nelem)    ); paw%vlocae=0.0d0
    allocate( paw%orb_2_l(norb,nelem)  ); paw%orb_2_l=-1000
    allocate( paw%orb_2_n(norb,nelem)  ); paw%orb_2_n=0

    call find_tag( g, "z_valence", cbuf )
    call get_dnum_from_string( cbuf, "z_valence=", pp%zion )
    pp%zps(ik)=nint(pp%zion)
    write(*,*) cbuf(1:20)

    call find_tag( g, "l_max", cbuf )
    call get_num_from_string( cbuf, "l_max=", pp%mlps(ik) )
    write(*,*) cbuf(1:20)

    allocate( paw%DIJ(norb,norb,nelem) ); paw%DIJ=0.0d0
    allocate( paw%QIJL(nrr,norb,norb,0:4,nelem) ); paw%QIJL=0.0d0
    allocate( paw%Q(norb*norb,nelem) ); paw%Q=0.0d0
    allocate( paw%MULTIPOLES(norb*norb*5,nelem) ); paw%MULTIPOLES=0

    call find_tag( g, "mesh_size", cbuf )
    call get_num_from_string( cbuf, "mesh_size=", nrr )
    write(*,*) cbuf(1:20)

    call find_tag( g, "number_of_wfc", cbuf )
    call get_num_from_string( cbuf, "number_of_wfc=", nwfc )
    write(*,*) cbuf(1:20)

    call find_tag( g, "number_of_proj", cbuf )
    call get_num_from_string( cbuf, "number_of_proj=", norb )
    write(*,*) cbuf(1:20)

    call find_tag( g, "<PP_R ", cbuf )
    read(g,*) pp%rad(1:nrr,ik)
    write(*,*) cbuf(1:20)

    paw%rad(1:nrr,ik)=pp%rad(1:nrr,ik)

    call find_tag( g, "<PP_RAB ", cbuf )
    read(g,*) paw%rab(1:nrr,ik)
    write(*,*) cbuf(1:20)

    call find_tag( g, "<PP_NLCC ", cbuf )
    read(g,*) paw%rhocps(1:nrr,ik)
    write(*,*) cbuf(1:20)

    call find_tag( g, "<PP_LOCAL ", cbuf )
    read(g,*) paw%vlocps(1:nrr,ik)
    write(*,*) cbuf(1:20)

    call find_tag( g, "<PP_NONLOCAL>", cbuf )
    write(*,*) cbuf(1:20)
    do i=1,norb
       read(g,'(a)') cbuf
       write(*,*) cbuf
       call get_num_from_string( cbuf, "angular_momentum=", l )
       paw%orb_2_l(i,ik)=l
       paw%orb_2_n(i,ik)=count(paw%orb_2_l(:,ik)==l)
       write(*,*) i,paw%orb_2_l(i,ik),paw%orb_2_n(i,ik)
       do
          read(g,'(a)') cbuf
          if ( index(cbuf,">") /= 0 ) exit
       end do
       read(g,*) paw%beta(1:nrr,i,ik)
       read(g,'(a)') cbuf
    end do ! i

    call find_tag( g, "<PP_DIJ ", cbuf ); write(*,*) cbuf(1:20)
    read(g,*) paw%DIJ(1:norb,1:norb,ik)

    call find_tag( g, "<PP_Q ", cbuf ); write(*,*) cbuf(1:20)
    read(g,*) paw%Q(:,ik)

    call find_tag( g, "<PP_MULTIPOLES ", cbuf ); write(*,*) cbuf(1:20)
    read(g,*) paw%MULTIPOLES(1:norb*norb*(pp%mlps(ik)+1),ik)

    do
       call find_tag( g, "<PP_QIJL", cbuf ); write(*,*) cbuf(1:20)
       call get_num_from_string( cbuf, "first_index=", i )
       call get_num_from_string( cbuf, "second_index=", j )
       call get_num_from_string( cbuf, "angular_momentum=", l )
       read(g,*) paw%QIJL(1:nrr,i,j,l,ik)
       read(g,'(a)') cbuf ;write(*,*) "aaaa",cbuf(1:20)
       read(g,'(a)') cbuf ;write(*,*) "aaa",cbuf(1:20)
       if ( index(cbuf,"</PP_AUGMENTATION") /= 0 ) then
          exit
       else
          backspace(g)
       end if
    end do

    call find_tag( g, "<PP_PSWFC", cbuf ); write(*,*) cbuf(1:20)
    do i=1,nwfc
       call find_tag( g, "<PP_CHI", cbuf ); write(*,*) cbuf(1:20)
       do
          read(g,'(a)') cbuf
          if ( index(cbuf,">") /= 0 ) exit
       end do
       read(g,*) pp%upp_f(1:nrr,i,ik)
       call find_tag( g, "</PP_CHI", cbuf )
    end do

    call find_tag( g, "<PP_FULL_WFC ", cbuf ); write(*,*) cbuf(1:20)
    do i=1,norb
       call find_tag( g, "<PP_AEWFC", cbuf ); write(*,*) cbuf(1:20)
       read(g,*) paw%aewf(1:nrr,i,ik)
       call find_tag( g, "</PP_AEWFC", cbuf )
    end do
    do i=1,norb
       call find_tag( g, "<PP_PSWFC", cbuf ); write(*,*) cbuf(1:20)
       read(g,*) paw%pswf(1:nrr,i,ik)
       call find_tag( g, "</PP_PSWFC", cbuf )
    end do

    call find_tag( g, "<PP_RHOATOM ", cbuf ); write(*,*) cbuf(1:20)
    read(g,*) paw%rhops(1:nrr,ik)

    call find_tag( g, "<PP_OCCUPATIONS ", cbuf ); write(*,*) cbuf(1:20)
    read(g,*) paw%occ(1:norb,ik)

    call find_tag( g, "<PP_AE_NLCC ", cbuf ); write(*,*) cbuf(1:20)
    read(g,*) paw%rhocae(1:nrr,ik)

    call find_tag( g, "<PP_AE_VLOC ", cbuf ); write(*,*) cbuf(1:20)
    read(g,*) paw%vlocae(1:nrr,ik)

!---

    pp%lref(ik) = ubound(pp%vpp,2)

    if ( paw%rad(1,ik) == 0.0d0 ) then
      stop 'stop@read_paw_upf(1)'
    else
      pp%mr(ik) = nrr + 1
      do i=2,nrr+1
        pp%rad(i,ik) = paw%rad(i-1,ik)
      end do
      pp%rad(1,ik) = 0.0d0
      do j=1,norb
        do i=1,nrr
          pp%vpp(i,j-1) = sqrt(0.5d0)*paw%beta(i,j,ik)
        end do
        pp%vpp(0,j-1) = pp%vpp(1,j-1)
      end do
      do i=1,nrr
        pp%vpp(i,pp%lref(ik)) = 0.5d0*paw%vlocps(i,ik)
      end do
      pp%vpp(0,pp%lref(ik)) = pp%vpp(1,pp%lref(ik))
    end if

    if ( all(paw%rab(2:nrr,ik)==paw%rab(1,ik)) ) then
      tmp = pp%rad(2,ik) - pp%rad(1,ik)
      do i = pp%mr(ik)+1, pp%nrmax
        pp%rad(i,ik) = (i-1)*tmp
      end do
    else
      tmp1 = log(paw%rad(1,ik))
      tmp2 = log(paw%rad(2,ik)) - log(paw%rad(1,ik))
      do i = pp%mr(ik)+1,pp%nrmax
        tmp = tmp1 + (i-2)*tmp2
        pp%rad(i,ik) = exp(tmp)
      end do
    end if

    pp%nrps(ik)=0
    do j=1,norb
      do i=nrr,1,-1
        if ( paw%beta(i,j,ik) /= 0.0d0 ) then
          pp%nrps(ik) = max( pp%nrps(ik), i )
          exit
        end if
      end do
    end do

    pp%rps(ik) = pp%rad(pp%nrps(ik),ik)
    write(*,*) 'pp%nrps,pp%rps',pp%nrps(ik), pp%rps(ik)

    write(*,*) "norb=",norb

    pp%nproj(:,ik)=0
    do j=1,norb
      l=paw%orb_2_l(j,ik)
      pp%nproj(l,ik) = pp%nproj(l,ik) + 1
    end do

    pp%mlps(ik) = maxval( paw%orb_2_l(:,ik) )

    write(*,*) "pp%mlps=",pp%mlps(ik)

    i=0
    l0=0
    do ll = 0, pp%mlps(ik)
    do l = l0, l0+pp%nproj(ll,ik)-1
      i=i+1
      write(*,*) i,ll
    end do
    end do

!---

    l=maxval( paw%orb_2_l(:,ik) )
    n=maxval( paw%orb_2_n(:,ik) )
    allocate( paw%nl_2_orb(0:l,n,ik) ); paw%nl_2_orb=0
    do i=1,norb
       l=paw%orb_2_l(i,ik)
       n=paw%orb_2_n(i,ik)
       paw%nl_2_orb(l,n,ik)=i
    end do

!---

    paw%rhoae(:,ik)=0.0d0
    do j=1,norb
       paw%rhoae(:,ik)=paw%rhoae(:,ik)+paw%occ(j,ik)*paw%aewf(:,j,ik)**2
    end do

!---

    allocate( work(nrr) ); work=0.0d0

    rewind 10
    do i=1,nrr
       write(10,'(1x,5g20.10)') pp%rad(i,ik),(paw%beta(i,j,ik),j=1,norb)
    end do

    rewind 11
    do i=1,nrr
       write(11,'(1x,5g20.10)') pp%rad(i,ik),paw%vlocps(i,ik),paw%vlocae(i,ik)
    end do

    rewind 12
    do i=1,nrr
       write(12,'(1x,5g20.10)') pp%rad(i,ik),(paw%pswf(i,j,ik),j=1,norb)
    end do

    rewind 13
    do i=1,nrr
       write(13,'(1x,5g20.10)') pp%rad(i,ik),(paw%aewf(i,j,ik),j=1,norb)
    end do

    rewind 14
    do i=1,nrr
       write(14,'(1x,5g20.10)') pp%rad(i,ik),paw%rhops(i,ik),paw%rhoae(i,ik)
    end do

    work=paw%aewf(:,1,ik)**2-paw%pswf(:,1,ik)**2
    rewind 15
    do i=1,nrr
       write(15,'(1x,5g20.10)') pp%rad(i,ik),qij(i),work(i)
    end do

    pi4=4.0d0*acos(-1.0d0)
    do l=0,2*pp%mlps(ik)
    do j=1,norb
    do i=1,norb
    write(*,'(1x,"qijl",3i3,es20.10)') i,j,l,sum(paw%QIJL(:,i,j,l,ik)*paw%rad(:,ik)**l*paw%rab(:,ik))
    end do
    end do
    end do
    write(*,*) '|aewf|**2',sum(paw%aewf(:,1,ik)**2*paw%rab(:,ik))
    write(*,*) '|aewf|**2',sum(paw%aewf(:,3,ik)**2*paw%rab(:,ik))
    write(*,*) '|pswf|**2',sum(paw%pswf(:,3,ik)**2*paw%rab(:,ik))
    write(*,*) "rhops",sum(paw%rhops(:,ik)*paw%rab(:,ik))
    write(*,*) "rhoae",sum(paw%rhoae(:,ik)*paw%rab(:,ik))
    write(*,*) "rhocps",sum(paw%rhocps(:,ik)*paw%rab(:,ik)) &
         ,sum(paw%rhocps(:,ik)*pp%rad(:,ik)**2*paw%rab(:,ik))*pi4
    write(*,*) "rhocae",sum(paw%rhocae(:,ik)*paw%rab(:,ik)) &
         ,sum(paw%rhocae(:,ik)*pp%rad(:,ik)**2*paw%rab(:,ik))*pi4

    do j=1,norb
       work=work+paw%occ(j,ik)*paw%pswf(:,j,ik)**2
    end do
    write(*,*) 'sum(density_from_pswf)',sum(work*paw%rab(:,ik))

    do j=1,norb-1
       write(*,*) '|beta*pswf|',sum(paw%beta(:,j+1,ik)*paw%pswf(:,j,ik)*paw%rab(:,ik))
    end do

    allocate( vh(nrr)  ); vh=0.0d0
    allocate( vxc(nrr) ); vxc=0.0d0

    call calc_hartree_radial_1d &
         ( pp%rad(1:nrr,ik),paw%rab(1:nrr,ik),paw%rhocae(1:nrr,ik),vh,tmp )

    exc=0.0d0
    do i=1,nrr
!       trho=paw%rhops(i,ik)/pp%rad(i,ik)**2/pi4 +paw%rhocps(i,ik)
       trho=paw%rhocae(i,ik)
       call PZxc( trho, tmp, tmp2 )
       exc=exc+tmp*trho*paw%rab(i,ik)*pp%rad(i,ik)**2*pi4
       vxc(i)=tmp+trho*tmp2
    end do
    write(*,*) "exc=",exc

    rewind 40
    do i=1,nrr
       write(40,'(1x,5g20.10)') pp%rad(i,ik),vh(i),vxc(i),paw%vlocps(i,ik)
    end do
    rewind 41
    do i=1,nrr
       write(41,'(1x,5g20.10)') pp%rad(i,ik),vh(i),vxc(i) &
            ,-6.0d0/pp%rad(i,ik),paw%vlocae(i,ik)
    end do

  end subroutine read_paw_upf

  subroutine get_num_from_string( buf, str, n )
    implicit none
    character(*),intent(in) :: buf, str
    integer,intent(out) :: n
    integer :: i, j, l, len_str
    character(8) :: str_out
    l = index( buf, str )
    if ( l /= 0 ) then
       len_str = len_trim( adjustl(str) )
       i=l+len_str+1
       j=index( buf(i:), '"' ) + i - 2
       str_out = buf(i:j)
       read(str_out,*) n
    end if
  end subroutine get_num_from_string

  subroutine get_dnum_from_string( buf, str, d )
    implicit none
    character(*),intent(in) :: buf, str
    real(8),intent(out) :: d
    integer :: i, j, l, len_str
    character(8) :: str_out
    l = index( buf, str )
    if ( l /= 0 ) then
       len_str = len_trim( adjustl(str) )
       i=l+len_str+1
       j=index( buf(i:), '"' ) + i - 2
       str_out = buf(i:j)
       read(str_out,*) d
    end if
  end subroutine get_dnum_from_string

  subroutine get_mr_norb_upf( g, mr, norb )
    implicit none
    integer,intent(in) :: g
    integer,intent(out) :: mr, norb
    character(200) :: cbuf
    do
       call find_tag( g, "mesh_size", cbuf )
       call get_num_from_string( cbuf, "mesh_size=", mr )
       call find_tag( g, "number_of_proj", cbuf )
       call get_num_from_string( cbuf, "number_of_proj=", norb )
       exit
    end do
  end subroutine get_mr_norb_upf

  subroutine find_tag( g, indx, cbuf )
    implicit none
    integer,intent(in) :: g
    character(*),intent(in) :: indx
    character(*),intent(out) :: cbuf
    integer :: l
    l=len(indx)
    do
       read(g,'(a)',END=99) cbuf
       cbuf=adjustl(cbuf)
       if ( cbuf(1:l) == indx ) return
    end do
99  continue
  end subroutine find_tag

end module read_paw_upf_module


text:
! cloning on any constant...
      program c3
      integer i, n, m
      parameter (n=5, m = n + 1 - 2 + 1)
      i = 4
      i = i + 1
      call clonee(i)
      call clonee(n)
      call clonee(m)
      call clonee(5)
      call clonee(5 + n - m)
      end

      subroutine clonee(i)
      integer i
      if (i.eq.5) print *, 'ok'
      end


text:
Real Function erhon(em1, em2, lb1, lb2, srt)
  Dimension arrayj(19), arrayl(19), arraym(19), arrayw(19), arrayb(19)
  Save
  Data arrayj/0.5, 1.5, 0.5, 0.5, 2.5, 2.5, 1.5, 0.5, 1.5, 3.5, 1.5, 0.5, 1.5, 0.5, 2.5, 0.5, 1.5, 2.5, 3.5/
  Data arrayl/1, 2, 0, 0, 2, 3, 2, 1, 1, 3, 1, 0, 2, 0, 3, 1, 1, 2, 3/
  Data arraym/1.44, 1.52, 1.535, 1.65, 1.675, 1.68, 1.70, 1.71, 1.72, 1.99, 1.60, 1.62, 1.70, 1.90, 1.905, 1.910, 1.86, 1.93, 1.95/
  Data arrayw/0.2, 0.125, 0.15, 0.15, 0.155, 0.125, 0.1, 0.11, 0.2, 0.29, 0.25, 0.16, 0.28, 0.15, 0.3, 0.22, 0.25, 0.25, 0.24/
  Data arrayb/0.15, 0.20, 0.05, 0.175, 0.025, 0.125, 0.1, 0.20, 0.53, 0.34, 0.05, 0.07, 0.15, 0.45, 0.45, 0.058, 0.08, 0.12, 0.08/
  pi = 3.1415926
  xs = 0
  Do ir = 1, 19
    If (ir<=8) Then
      If (((lb1*lb2==27 .And. (lb1==1 .Or. lb2==1)) .Or. (lb1*lb2==25*2 .And. (lb1==2 .Or. lb2==2))) .Or. ((lb1*lb2==-25 .And. (lb1==-1 .Or. lb2==-1)) .Or. (lb1*lb2==-27*2 .And. (lb1==-2 .Or. lb2==-2)))) branch = 0.
      If ((iabs(lb1*lb2)==26 .And. (iabs(lb1)==1 .Or. iabs(lb2)==1)) .Or. (iabs(lb1*lb2)==26*2 .And. (iabs(lb1)==2 .Or. iabs(lb2)==2))) branch = 1./3.
      If (((lb1*lb2==27*2 .And. (lb1==2 .Or. lb2==2)) .Or. (lb1*lb2==25 .And. (lb1==1 .Or. lb2==1))) .Or. ((lb1*lb2==-25*2 .And. (lb1==-2 .Or. lb2==-2)) .Or. (lb1*lb2==-27 .And. (lb1==-1 .Or. lb2==-1)))) branch = 2./3.
    Else
      If (((lb1*lb2==27 .And. (lb1==1 .Or. lb2==1)) .Or. (lb1*lb2==25*2 .And. (lb1==2 .Or. lb2==2))) .Or. ((lb1*lb2==-25 .And. (lb1==-1 .Or. lb2==-1)) .Or. (lb1*lb2==-27*2 .And. (lb1==-2 .Or. lb2==-2)))) branch = 1.
      If ((iabs(lb1*lb2)==26 .And. (iabs(lb1)==1 .Or. iabs(lb2)==1)) .Or. (iabs(lb1*lb2)==26*2 .And. (iabs(lb1)==2 .Or. iabs(lb2)==2))) branch = 2./3.
      If (((lb1*lb2==27*2 .And. (lb1==2 .Or. lb2==2)) .Or. (lb1*lb2==25 .And. (lb1==1 .Or. lb2==1))) .Or. ((lb1*lb2==-25*2 .And. (lb1==-2 .Or. lb2==-2)) .Or. (lb1*lb2==-27 .And. (lb1==-1 .Or. lb2==-1)))) branch = 1./3.
    End If
    xs0 = fdr(arraym(ir), arrayj(ir), arrayl(ir), arrayw(ir), arrayb(ir), srt, em1, em2)
    xs = xs + 1.3*pi*branch*xs0*(0.1973)**2
  End Do
  erhon = xs
  Return
End Function erhon

max_stars_repo_name: xiaohaijin/AMPT

max_issues_repo_name: xiaohaijin/AMPT

max_forks_repo_name: xiaohaijin/AMPT


text:
module tdd_jams
  implicit none
  private

  public :: say_hello
contains
  subroutine say_hello
    print *, "Hello, tdd_jams!"
  end subroutine say_hello
end module tdd_jams

max_stars_repo_name: sourceryinstitute/miniFAVOR

max_issues_repo_name: sourceryinstitute/miniFAVOR

max_forks_repo_name: sourceryinstitute/miniFAVOR


text:
MODULE m_fft3dxc
CONTAINS
   SUBROUTINE fft3dxc( &
      afft, bfft, fg3, &
      k1d, k2d, k3d, ng3, kimax, isn, &
      igfft1, igfft2, pgfft, nstr)

!***********************************************************
!                                                          *
! interface for fg3(star) -- fft --> (a,b)fft (r) (isn=+1) *
!         or (a,b)fft (r) -- fft --> fg3(star)    (isn=-1) *
!                                                          *
! dimension of (a,b)fft is (k1d x k2d x k3d)               *
! afft and bfft contain the real/imaginary part of the fft *
! igfft1(i)  is the pointer from the g-sphere to stars     *
! igfft2(i)  is the pointer from the g-sphere to fft-grid  *
! pgfft(i)   contains the phases of the g-vectors of sph.  *
!                                                          *
!***********************************************************

      USE m_fft_interface
      IMPLICIT NONE

      INTEGER :: k1d, k2d, k3d, ng3, kimax, isn
      INTEGER :: igfft1(0:k1d*k2d*k3d-1), igfft2(0:k1d*k2d*k3d-1)
      INTEGER :: nstr(ng3)
      COMPLEX pgfft(0:k1d*k2d*k3d-1)
      REAL ::    afft(0:k1d*k2d*k3d-1), bfft(0:k1d*k2d*k3d-1)
      COMPLEX :: fg3(ng3)

      INTEGER :: i, ifftd
      REAL :: scale, zero
      COMPLEX :: ctmp

      LOGICAL :: forw
      INTEGER :: length_zfft(3)
      complex :: zfft(0:k1d*k2d*k3d-1)

      ifftd = k1d*k2d*k3d
      zero = 0.0

      IF (isn > 0) THEN

         !  ---> put stars onto the fft-grid

         afft = 0.0
         bfft = 0.0
         DO i = 0, kimax - 1
            ctmp = fg3(igfft1(i))*pgfft(i)
            afft(igfft2(i)) = real(ctmp)
            bfft(igfft2(i)) = aimag(ctmp)
         ENDDO
      ENDIF

!---> now do the fft (isn=+1 : g -> r ; isn=-1 : r -> g)

      zfft = cmplx(afft, bfft)
      if (isn == -1) then
         forw = .TRUE.
      else
         forw = .FALSE.
      end if
      length_zfft(1) = k1d
      length_zfft(2) = k2d
      length_zfft(3) = k3d
      call fft_interface(3, length_zfft, zfft, forw, igfft2(0:kimax-1))
      afft = real(zfft)
      bfft = aimag(zfft)

      IF (isn < 0) THEN

         !  ---> collect stars from the fft-grid

         DO i = 1, ng3
            fg3(i) = cmplx(0.0, 0.0)
         ENDDO
         scale = 1.0/ifftd
         DO i = 0, kimax - 1
            fg3(igfft1(i)) = fg3(igfft1(i)) + CONJG(pgfft(i))* &
                             zfft(igfft2(i))
         ENDDO
         DO i = 1, ng3
            fg3(i) = scale*fg3(i)/nstr(i)
         ENDDO
      ENDIF

   END SUBROUTINE fft3dxc
END MODULE m_fft3dxc

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text: #redirect Users/GreggAlexander

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
c This module/program is part of DaliLite (c) L. Holm 1999,2008
c
	program wolf
c
c	f77 wolf.f
c
c	use as: wolf_grid | wolf
c
c       reads WOLFDB
c
	implicit none
	include 'parsizes.for'
        integer maxlines
        parameter (maxlines=1000000)
c       largest protein: 1y1vA 1426 residues, 69 SSEs, 4041 frame-lines
	integer*2 box(-20:20,-20:20,-20:20),link_list(maxlines)
c       parsizes.for has maxseg=200
        integer q_aseg(maxlines),q_bseg(maxlines),q_cseg(maxlines)
        character*1 q_atype(maxlines),q_ctype(maxlines)
        real q_mid(3,maxlines),q_dir(3,maxlines)
        character*5 query_cd
c
        character*5 cd,old_cd
        real mid(3),dir(3)
        character*1 atype,ctype
        integer aseg,bseg,cseg,nseg
c
        integer i,n,old_aseg,old_bseg
        integer gx,gy,gz
        real scores(maxseg,maxseg),best_seen
        logical less
        integer dx,dy,dz,gx0,gy0,gz0
        real score,cosi,r2,distance_squared,cosine
c
        real bestscore,bestperframe_score(maxseg,maxseg)
        integer best_iseg,best_jseg
        integer bestperframe_ijseg(2,maxseg,maxseg)
	character*80 filename
c
c	binary database
c
	write(*,*) 'enter WOLF90 database'
	read(*,500) filename
500	format(a80)

c       initialise box
c
        call initgrid(box,20,20,20)
c
c       read query from fort.1 (binary)
c
        n=1
10      read(1, end=19) query_cd,q_aseg(n),q_bseg(n),q_cseg(n),
     $          q_atype(n),q_ctype(n),gx,gy,gz,
     $          (q_mid(i,n),i=1,3),(q_dir(i,n),i=1,3)

c       load linked list of line-numbers in box
        i=box(gx,gy,gz)
        link_list(n)=i
        box(gx,gy,gz)=n
        n=n+1
        goto 10 ! query input loop
19      continue
c
c       read targets stream from WOLFDB (binary)
c
        old_cd='?????'
        old_aseg=0
        old_bseg=0
        nseg=0
        open(90,file=filename,form='unformatted',
     $          status='old')
20      read(90,end=29) cd,aseg,bseg,cseg,atype,ctype,
     $          gx0,gy0,gz0,(mid(i),i=1,3),(dir(i),i=1,3)

	if(aseg.gt.nseg) nseg=aseg
        if(bseg.gt.nseg) nseg=bseg
        if(cseg.gt.nseg) nseg=cseg
c       new target
        if(cd.ne.old_cd) then
                if(old_cd.ne.'?????')
     $		call output(query_cd,old_cd,nseg,
     $          bestperframe_ijseg,
     $          bestperframe_score)
                call init_scores(bestperframe_score,nseg)
                call init_counts(bestperframe_ijseg,nseg)
                call init_scores(scores,nseg)
                old_aseg=0
                old_bseg=0
        end if
        old_cd=cd
c       new target.frame
        if(aseg.ne.old_aseg.and.bseg.ne.old_bseg) then
          if(old_aseg.gt.0.and.old_bseg.gt.0) then
            call save_frame(scores,nseg,best_iseg,best_jseg,
     $          bestscore)
            bestperframe_ijseg(1,old_aseg,old_bseg)=best_iseg
            bestperframe_ijseg(2,old_aseg,old_bseg)=best_jseg
            bestperframe_score(old_aseg,old_bseg)=bestscore
          end if
        end if
        old_aseg=aseg
        old_bseg=bseg
c       cumulate comparison scores per query.frame
c       find best one-to-one segment match of target.segment
        less=aseg.lt.bseg
        best_seen=0.0
        do dx=-2,2
                gx=gx0+dx
		if(gx.lt.-20.or.gx.gt.20) goto 31
                do dy=-2,2
                        gy=gy0+dy
			if(gy.lt.-20.or.gy.gt.20) goto 32
                        do dz=-2,2
                              gz=gz0+dz
	if(gz.lt.-20.or.gz.gt.20) goto 33
                              i=box(gx,gy,gz)
                      do while(i.gt.0)
        if(atype.ne.q_atype(i)) goto 39
        if(ctype.ne.q_ctype(i)) goto 39
        if(less.and.q_aseg(i).ge.q_cseg(i)) goto 39
        if(.not.less.and.q_aseg(i).le.q_cseg(i)) goto 39
        r2=distance_squared(mid(1),mid(2),mid(3),
     $          q_mid(1,i),q_mid(2,i),q_mid(3,i))
        if(r2.gt.16.0) goto 39
        cosi=cosine(dir(1),dir(2),dir(3),
     $          q_dir(1,i),q_dir(2,i),q_dir(3,i))
        if(cosi.lt.0.5) goto 39
        score=1.0-r2/16.0+(cosi-0.5)*2
        if(score.gt.best_seen) best_seen=score
39                                      i=link_list(i)
                              end do
33                      end do
32              end do
31      end do
c       only allow one-to-one segment matching of target segment
        scores(aseg,bseg)=scores(aseg,bseg)+best_seen
        goto 20 ! WOLFDB input loop
29      close(90)
        if(old_aseg.ne.0.and.old_bseg.ne.0) then
          if(old_aseg.gt.0.and.old_bseg.gt.0) then
            call save_frame(scores,nseg,best_iseg,best_jseg,
     $          bestscore)
            bestperframe_ijseg(1,old_aseg,old_bseg)=best_iseg
            bestperframe_ijseg(2,old_aseg,old_bseg)=best_jseg
            bestperframe_score(old_aseg,old_bseg)=bestscore
          end if
        end if
        if(old_cd.ne.'?????') call output(query_cd,old_cd,nseg,
     $          bestperframe_ijseg,
     $          bestperframe_score)

550     format(8x,a5,3i5,2(1x,a1),3i5,6f10.1)

        end
c
c----------------------------------------------------------------------
c
        subroutine save_frame(scores,nseg,best_iseg,best_jseg,bestscore)
        implicit none
        include 'parsizes.for'
        real bestscore
        integer best_iseg,best_jseg
        real scores(maxseg,maxseg)
        integer nseg
c
        integer iseg,jseg

        bestscore=0.0
        best_iseg=0
        best_jseg=0
        do iseg=1,nseg
                do jseg=1,nseg
                        if(scores(iseg,jseg).gt.bestscore) then
                               bestscore=scores(iseg,jseg)
                               best_iseg=iseg
                               best_jseg=jseg
                        end if
                end do
        end do

        return
        end
c
c----------------------------------------------------------------------
c
        subroutine output(query_cd,target_cd,nseg,bestperframe_ijseg,
     $          bestperframe_score)
        implicit none
        include 'parsizes.for'
        character*5 query_cd,target_cd
        integer nseg,iseg,jseg
        integer bestperframe_ijseg(2,maxseg,maxseg)
        real bestperframe_score(maxseg,maxseg)
c
        real x,protcount
        integer lseg,kseg,useg,vseg,bestpair(4)
c
c       report best match per target_cd
c
	x=0.0
	useg=0
	vseg=0
        iseg=0
        jseg=0
	do lseg=1,nseg
	  do kseg=1,nseg
		if(bestperframe_score(lseg,kseg).gt.x) then
			x=bestperframe_score(lseg,kseg)
			useg=lseg
			vseg=kseg
                        iseg=bestperframe_ijseg(1,lseg,kseg)
                        jseg=bestperframe_ijseg(2,lseg,kseg)
		end if
	  end do
	end do
        write(*,500) query_cd,iseg,jseg,
     $          target_cd,useg,vseg,x

500     format(a5,2i5,1x,a5,2i5,f10.1)

        return
        end
c
c----------------------------------------------------------------------
c
	function cosine(z1,z2,z3,x1,x2,x3)
	implicit none
	real cosine,z1,z2,z3,x1,x2,x3
c
	real lz,lx,norm
c
	lz=sqrt(z1*z1+z2*z2+z3*z3)
	lx=sqrt(x1*x1+x2*x2+x3*x3)
	norm=max(1e-6,lz*lx)
	cosine=(z1*x1+z2*x2+z3*x3)/norm
c
	return
	end
c
c----------------------------------------------------------------------
c
	function distance_squared(a1,a2,a3,b1,b2,b3)
	implicit none
	real distance_squared,a1,a2,a3,b1,b2,b3
c
	distance_squared=(a1-b1)*(a1-b1)+(a2-b2)*(a2-b2)+(a3-b3)*(a3-b3)
c
	return
	end
c
c----------------------------------------------------------------------
c
	subroutine initgrid(box,boxdim1,boxdim2,boxdim3)
	implicit none
	integer boxdim1,boxdim2,boxdim3
	integer*2 box(-boxdim1:boxdim1,-boxdim2:boxdim2,-boxdim3:boxdim3)
c
	integer i,j,k
c
	do k=-boxdim3,boxdim3
		do j=-boxdim2,boxdim2
			do i=-boxdim1,boxdim1
				box(i,j,k)=0
			end do
		end do
	end do

	return
	end
c
c------------------------------------------------------------------------------
c
        subroutine init_scores(scores,nseg)
        implicit none
        include 'parsizes.for'
        integer nseg
        real scores(maxseg,maxseg)
c
        integer i,j

        do i=1,nseg
                do j=1,nseg
                        scores(i,j)=0.0
                end do
        end do

        return
        end
c
c------------------------------------------------------------------------------
c
        subroutine init_counts(scores,nseg)
        implicit none
        include 'parsizes.for'
        integer nseg
        integer scores(2,maxseg,maxseg)
c
        integer i,j,k

        do i=1,nseg
                do j=1,nseg
                        do k=1,2
                                scores(k,i,j)=0
                        end do
                end do
        end do

        return
        end
c
c------------------------------------------------------------------------------
c

max_stars_repo_name: gkrash/lcspromatch

max_issues_repo_name: gkrash/lcspromatch

max_forks_repo_name: gkrash/lcspromatch


subroutine bkgvar_rewgt(sfvar,vpvar,tvar,psvar,mype)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    bkgvar_rewgt   add flow dependence to variances
!   prgmmr: kleist           org: np22                date: 2007-07-03
!
! abstract: perform background error variance reweighting based on
!           flow dependence
!
! program history log:
!   2007-07-03  kleist
!   2008-06-05  safford - rm unused uses
!   2008-09-05  lueken - merged ed's changes into q1fy09 code
!   2009-04-15  wgu - added fpsproj option
!   2009-04-21  wgu - bug fix in routine smooth2d
!   2010-03-31  treadon - replace specmod components with sp_a structure
!   2011-10-09  wgu - add fut2ps to project unbalanced temp to surface pressure in static B modeling
!   2012-06-25  parrish - reorganize subroutine smooth2d so use one call in place of 4 calls.
!   2013-10-19  todling - metguess now holds background; count to fcount (count is intrisic function);
!                         protect against calls from non-FGAT run
!   2014-12-03  derber - restructure to do difference before sub2grid to
!                        optimize code
!
!   input argument list:
!     sfvar     - stream function variance
!     vpvar     - unbalanced velocity potential variance
!     tvar      - unbalanced temperature variance
!     psvar     - unbalanced surface pressure variance
!     mype      - integer mpi task id
!
!   output argument list:
!     sfvar     - reweighted stream function variance
!     vpvar     - reweighted unbalanced velocity potential variance
!     tvar      - reweighted unbalanced temperature variance
!     psvar     - reweighted unbalanced surface pressure variance
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_kind,i_kind,r_quad
  use constants, only: one,zero,two,zero_quad,tiny_r_kind
  use gridmod, only: nlat,nlon,nsig,lat2,lon2
  use guess_grids, only: nfldsig
  use mpimod, only: npe,mpi_comm_world,ierror,mpi_sum,mpi_rtype,mpi_max
  use balmod, only: agvz,wgvz,bvz,pput
  use berror, only: bkgv_rewgtfct,bkgv_write,fpsproj,fut2ps
  use gsi_metguess_mod, only: gsi_metguess_bundle
  use gsi_bundlemod, only: gsi_bundlegetpointer
  use mpeu_util, only: die
  implicit none

! Declare passed variables
  real(r_kind),dimension(lat2,lon2,nsig),intent(inout) :: sfvar,vpvar,tvar
  real(r_kind),dimension(lat2,lon2)     ,intent(inout) :: psvar
  integer(i_kind)                       ,intent(in   ) :: mype

! Declare local variables
  character(len=*),parameter::myname='bkgvar_rewgt'
  real(r_kind),dimension(lat2,lon2,nsig):: bald,balt
  real(r_kind),dimension(lat2,lon2,nsig):: delpsi,delchi,deltv
  real(r_kind),dimension(lat2,lon2):: delps,psresc,balps

  real(r_quad),dimension(nsig):: mean_dz,mean_dd,mean_dt
  real(r_quad) mean_dps,fcount

  real(r_kind),dimension(nsig,2,npe):: mean_dz0,mean_dd0,mean_dz1,mean_dd1,&
       mean_dt0,mean_dt1
  real(r_kind),dimension(2,npe):: mean_dps0,mean_dps1

  real(r_kind),dimension(nsig):: mean_dzout,mean_ddout,mean_dtout
  real(r_kind) mean_dpsout

  real(r_kind),dimension(nsig):: max_dz,max_dd,max_dt,max_dz0,max_dd0,&
       max_dt0,rmax_dz0,rmax_dd0,rmax_dt0
  real(r_kind) max_dps,max_dps0,rmax_dps0
  integer(i_kind) i,j,k,l,nsmth,mm1,nf,ier,istatus

  real(r_kind),dimension(:,:  ),pointer::ges_ps_01 =>NULL()
  real(r_kind),dimension(:,:  ),pointer::ges_ps_nf =>NULL()
  real(r_kind),dimension(:,:,:),pointer::ges_div_01=>NULL()
  real(r_kind),dimension(:,:,:),pointer::ges_div_nf=>NULL()
  real(r_kind),dimension(:,:,:),pointer::ges_vor_01=>NULL()
  real(r_kind),dimension(:,:,:),pointer::ges_vor_nf=>NULL()
  real(r_kind),dimension(:,:,:),pointer::ges_tv_01 =>NULL()
  real(r_kind),dimension(:,:,:),pointer::ges_tv_nf =>NULL()

! Initialize local arrays
  psresc=zero

  mean_dpsout=zero ; mean_dzout=zero ; mean_ddout=zero ; mean_dtout=zero
  mean_dps0  =zero ; mean_dz0  =zero ; mean_dd0  =zero ; mean_dt0  =zero
  mean_dps1  =zero ; mean_dz1  =zero ; mean_dd1  =zero ; mean_dt1  =zero

  max_dps=zero ; max_dps0=zero
  max_dz =zero ; max_dd  =zero ; max_dt=zero ; max_dz0=zero
  max_dd0=zero ; max_dt0 =zero
  balt   =zero ; bald    =zero ; balps =zero

! Set count to number of global grid points in quad precision
  fcount = float(nlat)*float(nlon)

! Set parameter for communication
  mm1=mype+1

! NOTES:
! This subroutine will only work if FGAT (more than one guess time level)
! is used.  For the current operational GDAS with a 6 hour time window,
! the default is to use sigf09-sigf03 to compute the delta variables
!
! Because of the FGAT component and vorticity/divergence issues, this is
! currently set up for global only
!
! No reweighting of ozone, cloud water, moisture yet
  if(nfldsig==1) return  ! not FGAT, nothing to do

! Get all pointers from met-guess
  nf=nfldsig
  ier=0
  call gsi_bundlegetpointer (gsi_metguess_bundle( 1),'ps'  ,ges_ps_01,    istatus)
  ier=ier+istatus
  call gsi_bundlegetpointer (gsi_metguess_bundle(nf),'ps'  ,ges_ps_nf,    istatus)
  ier=ier+istatus
  call gsi_bundlegetpointer (gsi_metguess_bundle( 1),'div' ,ges_div_01,   istatus)
  ier=ier+istatus
  call gsi_bundlegetpointer (gsi_metguess_bundle(nf),'div' ,ges_div_nf,   istatus)
  ier=ier+istatus
  call gsi_bundlegetpointer (gsi_metguess_bundle( 1),'vor' ,ges_vor_01,   istatus)
  ier=ier+istatus
  call gsi_bundlegetpointer (gsi_metguess_bundle(nf),'vor' ,ges_vor_nf,   istatus)
  ier=ier+istatus
  call gsi_bundlegetpointer (gsi_metguess_bundle( 1),'tv'  ,ges_tv_01,    istatus)
  ier=ier+istatus
  call gsi_bundlegetpointer (gsi_metguess_bundle(nf),'tv'  ,ges_tv_nf,    istatus)
  ier=ier+istatus
  if(ier/=0) call die(myname,'missing fields, ier= ', ier)

! Get stream function and velocity potential from guess vorticity and divergence
  call getpsichi(ges_vor_01,ges_vor_nf,delpsi)
  call getpsichi(ges_div_01,ges_div_nf,delchi)

! Get delta variables
  do k=1,nsig
     do j=1,lon2
        do i=1,lat2
           deltv(i,j,k) =ges_tv_nf(i,j,k)-ges_tv_01(i,j,k)
        end do
     end do
  end do
  do j=1,lon2
     do i=1,lat2
        delps(i,j) = ges_ps_nf(i,j)-ges_ps_01(i,j)
     end do
  end do

! Balanced surface pressure and velocity potential from delta stream function
  do k=1,nsig
     do j=1,lon2
        do i=1,lat2
           bald(i,j,k)=bvz(i,k)*delpsi(i,j,k)
        end do
     end do
  end do

! Balanced temperature from delta stream function
  do k=1,nsig
     do l=1,nsig
        do j=1,lon2
           do i=1,lat2
              balt(i,j,l)=balt(i,j,l)+agvz(i,l,k)*delpsi(i,j,k)
           end do
        end do
     end do
  end do

! Subtract off balanced parts
  do k=1,nsig
     do j=1,lon2
        do i=1,lat2
           deltv (i,j,k) = deltv (i,j,k) - balt(i,j,k)
           delchi(i,j,k) = delchi(i,j,k) - bald(i,j,k)
        end do
     end do
  end do

  if(fpsproj)then
     do k=1,nsig
        do j=1,lon2
           do i=1,lat2
              balps(i,j)=balps(i,j)+wgvz(i,k)*delpsi(i,j,k)
           end do
        end do
     end do
     if(fut2ps)then
        do k=1,nsig
           do j=1,lon2
              do i=1,lat2
                 balps(i,j)=balps(i,j)+pput(i,k)*deltv(i,j,k)
              end do
           end do
        end do
     endif
  else
     do j=1,lon2
        do i=1,lat2
           do k=1,nsig-1
              balps(i,j)=balps(i,j)+wgvz(i,k)*delpsi(i,j,k)
           end do
           balps(i,j)=balps(i,j)+wgvz(i,nsig)*delchi(i,j,1)
        end do
     end do
  endif

  do j=1,lon2
     do i=1,lat2
        delps(i,j) = delps(i,j) - balps(i,j)
     end do
  end do

! Convert to root mean square
  do k=1,nsig
     do j=1,lon2
        do i=1,lat2
           delpsi(i,j,k)=sqrt( delpsi(i,j,k)**two )
           delchi(i,j,k)=sqrt( delchi(i,j,k)**two )
           deltv (i,j,k)=sqrt( deltv (i,j,k)**two )
        end do
     end do
  end do
  do j=1,lon2
     do i=1,lat2
        delps(i,j)=sqrt( delps(i,j)**two )
     end do
  end do

! Smooth the delta fields before computing variance reweighting
  nsmth=8
  call smooth2d(delpsi,delchi,deltv,delps,nsig,nsmth,mype)

! Get global maximum and mean of each of the delta fields; while accounting for
! reproducibility of this kind of calculation
  mean_dz =zero_quad ; mean_dd=zero_quad ; mean_dt=zero_quad
  mean_dps=zero_quad

  do k=1,nsig
     do j=2,lon2-1
        do i=2,lat2-1
           mean_dz(k) = mean_dz(k) + delpsi(i,j,k)
           mean_dd(k) = mean_dd(k) + delchi(i,j,k)
           mean_dt(k) = mean_dt(k) + deltv (i,j,k)

           max_dz(k)=max(max_dz(k),delpsi(i,j,k))
           max_dd(k)=max(max_dd(k),delchi(i,j,k))
           max_dt(k)=max(max_dt(k),deltv (i,j,k))
        end do
     end do
  end do
  do j=2,lon2-1
     do i=2,lat2-1
        mean_dps = mean_dps + delps(i,j)
        max_dps  = max(max_dps,delps(i,j))
     end do
  end do

! Break quadruple precision into two double precision arrays
  do k=1,nsig
     mean_dz0(k,1,mm1) = mean_dz(k)
     mean_dz0(k,2,mm1) = mean_dz(k) - mean_dz0(k,1,mm1)
     mean_dd0(k,1,mm1) = mean_dd(k)
     mean_dd0(k,2,mm1) = mean_dd(k) - mean_dd0(k,1,mm1)
     mean_dt0(k,1,mm1) = mean_dt(k)
     mean_dt0(k,2,mm1) = mean_dt(k) - mean_dt0(k,1,mm1)
  end do
  mean_dps0(1,mm1) = mean_dps
  mean_dps0(2,mm1) = mean_dps - mean_dps0(1,mm1)

! Get task specific max and mean to every task
  call mpi_allreduce(mean_dz0,mean_dz1,nsig*2*npe,mpi_rtype,mpi_sum,mpi_comm_world,ierror)
  call mpi_allreduce(mean_dd0,mean_dd1,nsig*2*npe,mpi_rtype,mpi_sum,mpi_comm_world,ierror)
  call mpi_allreduce(mean_dt0,mean_dt1,nsig*2*npe,mpi_rtype,mpi_sum,mpi_comm_world,ierror)
  call mpi_allreduce(mean_dps0,mean_dps1,2*npe,mpi_rtype,mpi_sum,mpi_comm_world,ierror)

  call mpi_allreduce(max_dz,max_dz0,nsig,mpi_rtype,mpi_max,mpi_comm_world,ierror)
  call mpi_allreduce(max_dd,max_dd0,nsig,mpi_rtype,mpi_max,mpi_comm_world,ierror)
  call mpi_allreduce(max_dt,max_dt0,nsig,mpi_rtype,mpi_max,mpi_comm_world,ierror)
  call mpi_allreduce(max_dps,max_dps0,1,mpi_rtype,mpi_max,mpi_comm_world,ierror)

! Reintegrate quad precision number and sum over all mpi tasks
  mean_dz =zero_quad ; mean_dd=zero_quad ; mean_dt=zero_quad
  mean_dps=zero_quad
  do i=1,npe
     do k=1,nsig
        mean_dz(k) = mean_dz(k) + mean_dz1(k,1,i) + mean_dz1(k,2,i)
        mean_dd(k) = mean_dd(k) + mean_dd1(k,1,i) + mean_dd1(k,2,i)
        mean_dt(k) = mean_dt(k) + mean_dt1(k,1,i) + mean_dt1(k,2,i)
     end do
     mean_dps = mean_dps + mean_dps1(1,i) + mean_dps1(2,i)
  end do

! Divide by number of grid points to get the mean
  do k=1,nsig
     mean_dz(k)=mean_dz(k)/fcount
     mean_dd(k)=mean_dd(k)/fcount
     mean_dt(k)=mean_dt(k)/fcount
  end do
  mean_dps = mean_dps/fcount

! Load quad precision array back into double precision array for use
  do k=1,nsig
     mean_dzout(k)=mean_dz(k)
     mean_ddout(k)=mean_dd(k)
     mean_dtout(k)=mean_dt(k)
  end do
  mean_dpsout = mean_dps

! Take reciprocal of max values.  Check for tiny values
  do k=1,nsig
     rmax_dz0(k)=zero
     rmax_dd0(k)=zero
     rmax_dt0(k)=zero

     if (abs(max_dz0(k))>tiny_r_kind) rmax_dz0(k)=one/max_dz0(k)
     if (abs(max_dd0(k))>tiny_r_kind) rmax_dd0(k)=one/max_dd0(k)
     if (abs(max_dt0(k))>tiny_r_kind) rmax_dt0(k)=one/max_dt0(k)
  end do
  rmax_dps0=zero
  if (abs(max_dps0)>tiny_r_kind) rmax_dps0=one/max_dps0

! Get rescaling factor for each of the variables based on factor, mean, and max
!$omp parallel do schedule(dynamic,1) private(i,j,k)
  do k=1,nsig
     do j=1,lon2
        do i=1,lat2
           sfvar(i,j,k)=sfvar(i,j,k)* &
                (one + bkgv_rewgtfct*(delpsi(i,j,k)-mean_dzout(k))*rmax_dz0(k))
           vpvar(i,j,k)=vpvar(i,j,k)* &
                (one + bkgv_rewgtfct*(delchi(i,j,k)-mean_ddout(k))*rmax_dd0(k))
           tvar(i,j,k)=tvar(i,j,k)*  &
                (one + bkgv_rewgtfct*(deltv (i,j,k)-mean_dtout(k))*rmax_dt0(k))
        end do
     end do
  end do
  do j=1,lon2
     do i=1,lat2
        psvar(i,j)=psvar(i,j)* &
             (one + bkgv_rewgtfct*(delps(i,j)-mean_dpsout)*rmax_dps0)
     end do
  end do

! Smooth background error variances and write out grid
  nsmth=8
  call smooth2d(sfvar,vpvar,tvar,psvar,nsig,nsmth,mype)

  if (bkgv_write) call write_bkgvars_grid(sfvar,vpvar,tvar,psvar,mype)
  if(mype==0) write(6,*) 'bkgvar_rewgt: Flow-dependent feature on: nt=',nfldsig, ' minus nt= 1'

  return
end subroutine bkgvar_rewgt

subroutine getpsichi(vordiv1,vordiv2,dpsichi)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    getpsichi   compute psi and chi on subdomains
!   prgmmr: kleist           org: np22                date: 2007-07-03
!
! abstract: get stream function and velocity potential from vorticity
!           and divergence
!
! program history log:
!   2007-07-03  kleist
!   2008-06-05  safford - rm unused uses
!   2008-09-05  lueken - merged ed's changes into q1fy09 code
!   2010-04-01  treadon - move strip,reorder,reorder2 to gridmod
!   2012-06-25  parrish - replace strip,reorder,reorder2, mpi_all2allv with
!                         general_sub2grid/general_grid2sub.
!
!   input argument list:
!     vor       - vorticity on subdomains
!     div       - divergence on subdomains
!
!   output argument list:
!     psi       - stream function on subdomains
!     chi       - velocity potential on subdomains
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_kind,i_kind
  use constants, only: zero
  use gridmod, only: lat2,nsig,lon2,nlat,nlon,sp_a,grd_a
  use general_commvars_mod, only: g3
  use general_sub2grid_mod, only: general_sub2grid,general_grid2sub
  implicit none

! Declare passed variables
  real(r_kind),dimension(lat2,lon2,nsig),intent(in   ) :: vordiv1,vordiv2
  real(r_kind),dimension(lat2,lon2,nsig),intent(  out) :: dpsichi

! Declare local variables
  integer(i_kind) i,ii,j,k

  real(r_kind),dimension(lat2*lon2*nsig) :: vd1
  real(r_kind),dimension(g3%inner_vars,nlat,nlon,g3%kbegin_loc:g3%kend_alloc):: work1
  real(r_kind),dimension(sp_a%nc):: spc1

  ii=0
  do k=1,nsig
     do j=1,lon2
        do i=1,lat2
           ii=ii+1
           vd1(ii)=vordiv2(i,j,k)-vordiv1(i,j,k)
        end do
     end do
  end do
  call general_sub2grid(g3,vd1,work1)

! Perform scalar g2s on work array
!$omp parallel do schedule(dynamic,1) private(k,spc1)
  do k=g3%kbegin_loc,g3%kend_loc
     spc1=zero
     call general_g2s0(grd_a,sp_a,spc1,work1(1,:,:,k))

! Inverse laplacian
     do i=2,sp_a%ncd2
        spc1(2*i-1)=spc1(2*i-1)/(-sp_a%enn1(i))
        spc1(2*i)=spc1(2*i)/(-sp_a%enn1(i))
     end do
     spc1(1)=zero
     spc1(2)=zero

     work1(1,:,:,k)=zero
     call general_s2g0(grd_a,sp_a,spc1,work1(1,:,:,k))

  end do

! Get psi/chi back on subdomains
  call general_grid2sub(g3,work1,vd1)
  ii=0
  do k=1,nsig
     do j=1,lon2
        do i=1,lat2
           ii=ii+1
           dpsichi(i,j,k)=vd1(ii)
        end do
     end do
  end do

  return
end subroutine getpsichi

subroutine smooth2d(subd1,subd2,subd3,subd4,nlevs,nsmooth,mype)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    smooth2d    perform nine point smoother to interior
!   prgmmr: kleist           org: np22                date: 2007-07-03
!
! abstract: perform communication necessary, then apply simple nine
!           point smoother to interior of domain
!
! program history log:
!   2007-07-03  kleist
!   2012-06-25  parrish - create new version of smooth2d using general_sub2grid/general_grid2sub.
!   2012-06-25  parrish - remove following subroutine scatter_stuff2 (no longer needed).
!
!   input argument list:
!     subd      - field to be smoothed on subdomains
!     nlevs     - number of levels to perform smoothing
!     nsmooth   - number of times to perform smoother
!     mype      - mpi integer task ik
!
!   output argument list:
!     subd      - smoothed field
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use gridmod,only: nlat,nlon,lat2,lon2
  use kinds,only: r_kind,i_kind
  use constants, only: zero,half,one,two,three,four
  use general_commvars_mod, only: g33p1
  use general_sub2grid_mod, only: general_sub2grid,general_grid2sub
  implicit none

  integer(i_kind),intent(in   ) :: nlevs,mype,nsmooth
  real(r_kind)   ,intent(inout) :: subd1(lat2,lon2,nlevs)
  real(r_kind)   ,intent(inout) :: subd2(lat2,lon2,nlevs)
  real(r_kind)   ,intent(inout) :: subd3(lat2,lon2,nlevs)
  real(r_kind)   ,intent(inout) :: subd4(lat2,lon2)

  real(r_kind),dimension(g33p1%inner_vars,lat2,lon2,g33p1%num_fields)::worksub
  real(r_kind),dimension(g33p1%inner_vars,nlat,nlon,g33p1%kbegin_loc:g33p1%kend_alloc)::grd
  real(r_kind),dimension(nlat,0:nlon+1):: grd2
  real(r_kind) corn,cent,side,temp,c2,c3,c4,rterm,sums,sumn
  integer i,j,k,kk,n

! Get subdomains on the grid
  kk=0
  do k=1,nlevs
     kk=kk+1
     do j=1,lon2
        do i=1,lat2
           worksub(1,i,j,kk)=subd1(i,j,k)
        end do
     end do
  end do
  do k=1,nlevs
     kk=kk+1
     do j=1,lon2
        do i=1,lat2
           worksub(1,i,j,kk)=subd2(i,j,k)
        end do
     end do
  end do
  do k=1,nlevs
     kk=kk+1
     do j=1,lon2
        do i=1,lat2
           worksub(1,i,j,kk)=subd3(i,j,k)
        end do
     end do
  end do
  kk=kk+1
  if(kk/=g33p1%num_fields) then
     if(mype==0) write(6,*)' 3*nlevs+1 /= g33p1%num_fields in smooth2d, program stops'
     call stop2(99)
  end if
     do j=1,lon2
        do i=1,lat2
           worksub(1,i,j,kk)=subd4(i,j)
        end do
     end do

  call general_sub2grid(g33p1,worksub,grd)

! Set weights for nine point smoother
  corn=0.3_r_kind
  cent=one
  side=half
  c4=four
  c3=three
  c2=two
  rterm = one/(cent + c4*side + c4*corn)

  do k=g33p1%kbegin_loc,g33p1%kend_loc

! Do nsmooth number of passes over the 9pt smoother
        do n=1,nsmooth

! Load grd2 which is used in computing the smoothed fields
           do j=1,nlon
              do i=1,nlat
                 grd2(i,j)=grd(1,i,j,k)
              end do
           end do

! Load longitude wrapper rows
           do i=1,nlat
              grd2(i,0)=grd(1,i,nlon,k)
              grd2(i,nlon+1)=grd(1,i,1,k)
           end do

! special treatment for near-poles
           sumn=zero
           sums=zero
           do i=1,nlon
              sumn=sumn+grd2(nlat-1,i)
              sums=sums+grd2(2,i)
           end do
           sumn=sumn/(real(nlon,r_kind))
           sums=sums/(real(nlon,r_kind))
           do i=0,nlon+1
              grd2(nlat,i)=sumn
              grd2(1,i)=sums
           end do
! Perform smoother on interior
           do j=1,nlon
              do i=2,nlat-1
                 temp = cent*grd2(i,j) + side*(grd2(i+1,j) + &
                    grd2(i-1,j) + grd2(i,j+1) + grd2(i,j-1)) + &
                    corn*(grd2(i+1,j+1) + grd2(i+1,j-1) + grd2(i-1,j-1) + &
                    grd2(i-1,j+1))
                 grd(1,i,j,k) = temp*rterm
              end do
           end do
        end do    ! n smooth number of passes
  end do  ! k levs

! Get back on subdomains
  call general_grid2sub(g33p1,grd,worksub)

  kk=0
  do k=1,nlevs
     kk=kk+1
     do j=1,lon2
        do i=1,lat2
           subd1(i,j,k)=worksub(1,i,j,kk)
        end do
     end do
  end do
  do k=1,nlevs
     kk=kk+1
     do j=1,lon2
        do i=1,lat2
           subd2(i,j,k)=worksub(1,i,j,kk)
        end do
     end do
  end do
  do k=1,nlevs
     kk=kk+1
     do j=1,lon2
        do i=1,lat2
           subd3(i,j,k)=worksub(1,i,j,kk)
        end do
     end do
  end do
  kk=kk+1
  if(kk/=g33p1%num_fields) then
     if(mype==0) write(6,*)' 3*nlevs+1 /= g33p1%num_fields in smooth2d, program stops'
     call stop2(99)
  end if
     do j=1,lon2
        do i=1,lat2
           subd4(i,j)=worksub(1,i,j,kk)
        end do
     end do

  return
end subroutine smooth2d

subroutine gather_stuff2(f,g,mype,outpe)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    gather_stuff2    gather subdomains onto global slabs
!   prgmmr: kleist           org: np22                date: 2007-07-03
!
! abstract: perform communication necessary to gather subdomains to
!           global slabs
!
! program history log:
!   2007-07-03  kleist
!   2010-04-01  treadon - move strip and reorder to gridmod
!   2010-06-02  kokron - protect my reorder & copy to work on outpe only
!
!   input argument list:
!     f        - field on subdomains
!     mype     - mpi integer task id
!     outpe    - task to output global slab onto
!
!   output argument list:
!     g        - global slab on task outpe
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_kind,i_kind
  use mpimod, only: mpi_rtype,mpi_comm_world,ierror
  use general_commvars_mod, only: g1
  implicit none

  integer(i_kind),intent(in   ) :: mype,outpe
  real(r_kind)   ,intent(in   ) :: f(g1%lat2,g1%lon2)
  real(r_kind)   ,intent(  out) :: g(g1%nlat,g1%nlon)

  real(r_kind) fsm(g1%lat1,g1%lon1)
  real(r_kind),allocatable:: tempa(:)
  integer(i_kind) i,isize,j,ioffset,ilat,jlon,n,i0,j0,iloc,iskip

  isize=max(g1%iglobal,g1%itotsub)
  allocate(tempa(isize))

! Strip off endpoints of input array on subdomains

  do j=2,g1%lon2-1
     j0=j-1
     do i=2,g1%lat2-1
        i0=i-1
        fsm(i0,j0)=f(i,j)
     end do
  end do
  call mpi_gatherv(fsm,g1%ijn(mype+1),mpi_rtype, &
                  tempa,g1%ijn,g1%displs_g,mpi_rtype,outpe,mpi_comm_world,ierror)

  if(mype==outpe) then
     iskip=0
     iloc=0
     do n=1,g1%npe
        ioffset=iskip
        do i=1,g1%ijn(n)
           iloc=iloc+1
           ilat=g1%ltosi(iloc)
           jlon=g1%ltosj(iloc)
           g(ilat,jlon)=tempa(i+ioffset)
        end do
        iskip=iskip+g1%ijn(n)
     end do
  endif

  deallocate(tempa)

end subroutine gather_stuff2