text:
************************************************************************
*     This is the W- routine                                           *
************************************************************************
      subroutine qqb_wmbjet_gs(p,msq)
************************************************************************
*     Author: J.M. Campbell                                            *
*     January, 2004.                                                   *
************************************************************************
c---Matrix element SUBTRACTION squared averaged over initial colors and spins
c     q(-p1)+b(-p2) -->  W + f(p5) + b(p6) + g(p7)
c                        |
c                        --> nu(p3) + e^+(p4)
c
c--- all momenta are incoming
c--- Extended to include charm quark production via the variable "flav"

      implicit none
      include 'constants.f'
      include 'ptilde.f'
      include 'qqgg.f'
      include 'heavyflav.f'

      integer n16_2,n17_2,n15_2,n26_1,
     & n27_1,n25_1,n57_6,n67_5,
     & n17_5,n57_1,n17_6,n67_1,n27_6,n67_2,n27_5,n57_2
      parameter(
     & n16_2=1,n17_2=2,n15_2=3,n26_1=4,
     & n27_1=5,n25_1=6,n57_6=7,n67_5=8,
     & n17_5=9,n57_1=9,
     & n17_6=10,n67_1=10,
     & n27_6=11,n67_2=11,
     & n27_5=12,n57_2=12)

      integer j,k
c --- remember: nd will count the dipoles
      integer nd,isub
c--- slightly obtuse notation, to simplify declaration lines
      double precision p(mxpart,4),msq(maxd,fn:nf,fn:nf)
      double precision
     & msq17_2(fn:nf,fn:nf),msq27_1(fn:nf,fn:nf),
     & msq15_2(fn:nf,fn:nf),msq25_1(fn:nf,fn:nf),
     & msq16_2(fn:nf,fn:nf),msq26_1(fn:nf,fn:nf),
     & msq17_5(fn:nf,fn:nf),msq17_6(fn:nf,fn:nf),
     & msq27_5(fn:nf,fn:nf),msq27_6(fn:nf,fn:nf),
     & msq57_6(fn:nf,fn:nf),msq67_5(fn:nf,fn:nf),
     & dummy(fn:nf,fn:nf),dummyv(fn:nf,fn:nf),
     & sub17_2(4),sub27_1(4),sub15_2(4),sub25_1(4),
     & sub16_2(4),sub26_1(4),
     & sub17_5(4),sub57_1(4),sub27_5(4),sub57_2(4),
     & sub17_6(4),sub67_1(4),sub27_6(4),sub67_2(4),
     & sub57_6(4),sub67_5(4),dsubv

      double precision
     & m17_2(0:2,fn:nf,fn:nf),m27_1(0:2,fn:nf,fn:nf),
     & m16_2(0:2,fn:nf,fn:nf),m26_1(0:2,fn:nf,fn:nf),
     & m15_2(0:2,fn:nf,fn:nf),m25_1(0:2,fn:nf,fn:nf),
     & m57_6(0:2,fn:nf,fn:nf),m67_5(0:2,fn:nf,fn:nf),
     & m17_5(0:2,fn:nf,fn:nf),
     & m17_6(0:2,fn:nf,fn:nf),
     & m27_5(0:2,fn:nf,fn:nf),
     & m27_6(0:2,fn:nf,fn:nf),
     & mqq(0:2,fn:nf,fn:nf)

      integer,parameter::jj(-nf:nf)=(/-1,-2,-1,-2,-1,0,1,2,1,2,1/)
      integer,parameter::kk(-nf:nf)=(/-1,-2,-1,-2,-1,0,1,2,1,2,1/)

      common/isub/isub
      external qqb_wbjet,donothing_gvec

      if (isub .eq. 2) then
        ndmax=6
      else
        ndmax=12
      endif

c-- initialize the matrix elements to zero
      do j=-nf,nf
      do k=-nf,nf
      do nd=1,ndmax
        msq(nd,j,k)=0d0
      enddo
      enddo
      enddo

c--- calculate all the dipoles
c--- the dipole number relates the matrix elements to the transformed
c--- momenta, used in realint to perform cuts and clustering

c--- NEW ATTEMPT - mqq plays the role of a dummy variable here
      call dips(n16_2,p,1,6,2,sub16_2,dsubv,msq16_2,dummyv,
     . qqb_wbjet,donothing_gvec)
      call storedip_mass(m16_2,mqq)
      call dips(n17_2,p,1,7,2,sub17_2,dsubv,msq17_2,dummyv,
     . qqb_wbjet,donothing_gvec)
      call storedip_mass(m17_2,mqq)
      call dips(n15_2,p,1,5,2,sub15_2,dsubv,msq15_2,dummyv,
     . qqb_wbjet,donothing_gvec)
      call storedip_mass(m15_2,mqq)
      call dips(n26_1,p,2,6,1,sub26_1,dsubv,msq26_1,dummyv,
     . qqb_wbjet,donothing_gvec)
      call storedip_mass(m26_1,mqq)
      call dips(n27_1,p,2,7,1,sub27_1,dsubv,msq27_1,dummyv,
     . qqb_wbjet,donothing_gvec)
      call storedip_mass(m27_1,mqq)
      call dips(n25_1,p,2,5,1,sub25_1,dsubv,msq25_1,dummyv,
     . qqb_wbjet,donothing_gvec)
      call storedip_mass(m25_1,mqq)

      if (isub .eq. 1) then
      call dips(n57_6,p,5,7,6,sub57_6,dsubv,msq57_6,dummyv,
     . qqb_wbjet,donothing_gvec)
      call storedip_mass(m57_6,mqq)
      call dips(n67_5,p,6,7,5,sub67_5,dsubv,msq67_5,dummyv,
     . qqb_wbjet,donothing_gvec)
      call storedip_mass(m67_5,mqq)

      call dips(n17_5,p,1,7,5,sub17_5,dsubv,msq17_5,dummyv,
     . qqb_wbjet,donothing_gvec)
      call storedip_mass(m17_5,mqq)
      call dips(n57_1,p,5,7,1,sub57_1,dsubv,dummy,dummyv,
     . qqb_wbjet,donothing_gvec)

      call dips(n17_6,p,1,7,6,sub17_6,dsubv,msq17_6,dummyv,
     . qqb_wbjet,donothing_gvec)
      call storedip_mass(m17_6,mqq)
      call dips(n67_1,p,6,7,1,sub67_1,dsubv,dummy,dummyv,
     . qqb_wbjet,donothing_gvec)

      call dips(n27_6,p,2,7,6,sub27_6,dsubv,msq27_6,dummyv,
     . qqb_wbjet,donothing_gvec)
      call storedip_mass(m27_6,mqq)
      call dips(n67_2,p,6,7,2,sub67_2,dsubv,dummy,dummyv,
     . qqb_wbjet,donothing_gvec)

      call dips(n27_5,p,2,7,5,sub27_5,dsubv,msq27_5,dummyv,
     . qqb_wbjet,donothing_gvec)
      call storedip_mass(m27_5,mqq)
      call dips(n57_2,p,5,7,2,sub57_2,dsubv,dummy,dummyv,
     . qqb_wbjet,donothing_gvec)
      endif

c--- fill the dipole contributions
      do j=-nf,nf
      do k=-nf,nf

      if ((abs(j) .ne. flav) .and. (abs(k) .ne. flav)) goto 97
      if ((abs(j) .eq. flav) .and. (abs(k) .eq. flav)) goto 97
c--- so that either abs(j) or abs(k) = flav (but not both).

c--- contributions with 1 b-quark
      if (isub .eq. 1) then

c--- Note that each term must be examined individually to decide
c--- whether the leading colour term is proportional to (1) or (2)
      if ((j .gt. 0) .and. (k.gt.0)) then
      if (j .eq. +flav) then
c---QQ - process bu -> bd
      msq(n57_6,j,k)=msq(n57_6,j,k)
     . +sub57_6(qq)
     . *((xn+1d0/xn)*m57_6(0,j,k)+two*(m57_6(1,j,k)+m57_6(2,j,k))/xn)
      msq(n67_5,j,k)=msq(n67_5,j,k)
     . +sub67_5(qq)
     . *((xn+1d0/xn)*m67_5(0,j,k)+two*(m67_5(1,j,k)+m67_5(2,j,k))/xn)
      msq(n17_5,j,k)=msq(n17_5,j,k)
     . +sub17_5(qq)
     . *((xn-two/xn)*m17_5(2,j,k)-m17_5(0,j,k)/xn-m17_5(1,j,k)/xn)
     . +sub57_1(qq)
     . *((xn-two/xn)*m17_5(2,j,k)-m17_5(0,j,k)/xn-m17_5(1,j,k)/xn)
      msq(n17_6,j,k)=msq(n17_6,j,k)
     . +sub17_6(qq)
     . *((xn-two/xn)*m17_6(1,j,k)-m17_6(0,j,k)/xn-m17_6(2,j,k)/xn)
     . +sub67_1(qq)
     . *((xn-two/xn)*m17_6(1,j,k)-m17_6(0,j,k)/xn-m17_6(2,j,k)/xn)
      msq(n27_6,j,k)=msq(n27_6,j,k)
     . +sub27_6(qq)
     . *((xn-two/xn)*m27_6(2,j,k)-m27_6(0,j,k)/xn-m27_6(1,j,k)/xn)
     . +sub67_2(qq)
     . *((xn-two/xn)*m27_6(2,j,k)-m27_6(0,j,k)/xn-m27_6(1,j,k)/xn)
      msq(n27_5,j,k)=msq(n27_5,j,k)
     . +sub27_5(qq)
     . *((xn-two/xn)*m27_5(1,j,k)-m27_5(0,j,k)/xn-m27_5(2,j,k)/xn)
     . +sub57_2(qq)
     . *((xn-two/xn)*m27_5(1,j,k)-m27_5(0,j,k)/xn-m27_5(2,j,k)/xn)
      msq(n17_2,j,k)=msq(n17_2,j,k)
     . +sub17_2(qq)
     . *((xn+1d0/xn)*m17_2(0,j,k)+two*(m17_2(1,j,k)+m17_2(2,j,k))/xn)
      msq(n27_1,j,k)=msq(n27_1,j,k)
     . +sub27_1(qq)
     . *((xn+1d0/xn)*m27_1(0,j,k)+two*(m27_1(1,j,k)+m27_1(2,j,k))/xn)
      else
c---QQ - process ub -> bd (obtained from above by switching p1 and p2)
c--- corrected this one so far (for 1 and 2 colour structure switch)
      msq(n57_6,j,k)=msq(n57_6,j,k)
     . +sub57_6(qq)
     . *((xn+1d0/xn)*m57_6(0,j,k)+two*(m57_6(2,j,k)+m57_6(1,j,k))/xn)
      msq(n67_5,j,k)=msq(n67_5,j,k)
     . +sub67_5(qq)
     . *((xn+1d0/xn)*m67_5(0,j,k)+two*(m67_5(2,j,k)+m67_5(1,j,k))/xn)
      msq(n27_5,j,k)=msq(n27_5,j,k)
     . +sub27_5(qq)
     . *((xn-two/xn)*m27_5(1,j,k)-m27_5(0,j,k)/xn-m27_5(2,j,k)/xn)
     . +sub57_2(qq)
     . *((xn-two/xn)*m27_5(1,j,k)-m27_5(0,j,k)/xn-m27_5(2,j,k)/xn)
      msq(n27_6,j,k)=msq(n27_6,j,k)
     . +sub27_6(qq)
     . *((xn-two/xn)*m27_6(2,j,k)-m27_6(0,j,k)/xn-m27_6(1,j,k)/xn)
     . +sub67_2(qq)
     . *((xn-two/xn)*m27_6(2,j,k)-m27_6(0,j,k)/xn-m27_6(1,j,k)/xn)
      msq(n17_6,j,k)=msq(n17_6,j,k)
     . +sub17_6(qq)
     . *((xn-two/xn)*m17_6(1,j,k)-m17_6(0,j,k)/xn-m17_6(2,j,k)/xn)
     . +sub67_1(qq)
     . *((xn-two/xn)*m17_6(1,j,k)-m17_6(0,j,k)/xn-m17_6(2,j,k)/xn)
      msq(n17_5,j,k)=msq(n17_5,j,k)
     . +sub17_5(qq)
     . *((xn-two/xn)*m17_5(2,j,k)-m17_5(0,j,k)/xn-m17_5(1,j,k)/xn)
     . +sub57_1(qq)
     . *((xn-two/xn)*m17_5(2,j,k)-m17_5(0,j,k)/xn-m17_5(1,j,k)/xn)
      msq(n17_2,j,k)=msq(n17_2,j,k)
     . +sub17_2(qq)
     . *((xn+1d0/xn)*m17_2(0,j,k)+two*(m17_2(2,j,k)+m17_2(1,j,k))/xn)
      msq(n27_1,j,k)=msq(n27_1,j,k)
     . +sub27_1(qq)
     . *((xn+1d0/xn)*m27_1(0,j,k)+two*(m27_1(2,j,k)+m27_1(1,j,k))/xn)
      endif
      elseif ((j .lt. 0) .and. (k.lt.0)) then
      if (j .eq. -flav) then
c--- QbarQbar - process bd -> bu (obtained from bu -> bd by switching
c---                               p1 and p5, p2 and p6)
      msq(n57_6,j,k)=msq(n57_6,j,k)
     . +sub57_6(qq)
     . *((xn+1d0/xn)*m57_6(0,j,k)+two*(m57_6(1,j,k)+m57_6(2,j,k))/xn)
      msq(n67_5,j,k)=msq(n67_5,j,k)
     . +sub67_5(qq)
     . *((xn+1d0/xn)*m67_5(0,j,k)+two*(m67_5(1,j,k)+m67_5(2,j,k))/xn)
      msq(n17_5,j,k)=msq(n17_5,j,k)
     . +sub17_5(qq)
     . *((xn-two/xn)*m17_5(2,j,k)-m17_5(0,j,k)/xn-m17_5(1,j,k)/xn)
     . +sub57_1(qq)
     . *((xn-two/xn)*m17_5(2,j,k)-m17_5(0,j,k)/xn-m17_5(1,j,k)/xn)
      msq(n17_6,j,k)=msq(n17_6,j,k)
     . +sub17_6(qq)
     . *((xn-two/xn)*m17_6(1,j,k)-m17_6(0,j,k)/xn-m17_6(2,j,k)/xn)
     . +sub67_1(qq)
     . *((xn-two/xn)*m17_6(1,j,k)-m17_6(0,j,k)/xn-m17_6(2,j,k)/xn)
      msq(n27_6,j,k)=msq(n27_6,j,k)
     . +sub27_6(qq)
     . *((xn-two/xn)*m27_6(2,j,k)-m27_6(0,j,k)/xn-m27_6(1,j,k)/xn)
     . +sub67_2(qq)
     . *((xn-two/xn)*m27_6(2,j,k)-m27_6(0,j,k)/xn-m27_6(1,j,k)/xn)
      msq(n27_5,j,k)=msq(n27_5,j,k)
     . +sub27_5(qq)
     . *((xn-two/xn)*m27_5(1,j,k)-m27_5(0,j,k)/xn-m27_5(2,j,k)/xn)
     . +sub57_2(qq)
     . *((xn-two/xn)*m27_5(1,j,k)-m27_5(0,j,k)/xn-m27_5(2,j,k)/xn)
      msq(n17_2,j,k)=msq(n17_2,j,k)
     . +sub17_2(qq)
     . *((xn+1d0/xn)*m17_2(0,j,k)+two*(m17_2(1,j,k)+m17_2(2,j,k))/xn)
      msq(n27_1,j,k)=msq(n27_1,j,k)
     . +sub27_1(qq)
     . *((xn+1d0/xn)*m27_1(0,j,k)+two*(m27_1(1,j,k)+m27_1(2,j,k))/xn)
      else
c---QbarQbar - process db -> bu (obtained from above by switching p1 and p2)
      msq(n57_6,j,k)=msq(n57_6,j,k)
     . +sub57_6(qq)
     . *((xn+1d0/xn)*m57_6(0,j,k)+two*(m57_6(2,j,k)+m57_6(1,j,k))/xn)
      msq(n67_5,j,k)=msq(n67_5,j,k)
     . +sub67_5(qq)
     . *((xn+1d0/xn)*m67_5(0,j,k)+two*(m67_5(2,j,k)+m67_5(1,j,k))/xn)
      msq(n27_5,j,k)=msq(n27_5,j,k)
     . +sub27_5(qq)
     . *((xn-two/xn)*m27_5(1,j,k)-m27_5(0,j,k)/xn-m27_5(2,j,k)/xn)
     . +sub57_2(qq)
     . *((xn-two/xn)*m27_5(1,j,k)-m27_5(0,j,k)/xn-m27_5(2,j,k)/xn)
      msq(n27_6,j,k)=msq(n27_6,j,k)
     . +sub27_6(qq)
     . *((xn-two/xn)*m27_6(2,j,k)-m27_6(0,j,k)/xn-m27_6(1,j,k)/xn)
     . +sub67_2(qq)
     . *((xn-two/xn)*m27_6(2,j,k)-m27_6(0,j,k)/xn-m27_6(1,j,k)/xn)
      msq(n17_6,j,k)=msq(n17_6,j,k)
     . +sub17_6(qq)
     . *((xn-two/xn)*m17_6(1,j,k)-m17_6(0,j,k)/xn-m17_6(2,j,k)/xn)
     . +sub67_1(qq)
     . *((xn-two/xn)*m17_6(1,j,k)-m17_6(0,j,k)/xn-m17_6(2,j,k)/xn)
      msq(n17_5,j,k)=msq(n17_5,j,k)
     . +sub17_5(qq)
     . *((xn-two/xn)*m17_5(2,j,k)-m17_5(0,j,k)/xn-m17_5(1,j,k)/xn)
     . +sub57_1(qq)
     . *((xn-two/xn)*m17_5(2,j,k)-m17_5(0,j,k)/xn-m17_5(1,j,k)/xn)
      msq(n17_2,j,k)=msq(n17_2,j,k)
     . +sub17_2(qq)
     . *((xn+1d0/xn)*m17_2(0,j,k)+two*(m17_2(2,j,k)+m17_2(1,j,k))/xn)
      msq(n27_1,j,k)=msq(n27_1,j,k)
     . +sub27_1(qq)
     . *((xn+1d0/xn)*m27_1(0,j,k)+two*(m27_1(2,j,k)+m27_1(1,j,k))/xn)
      endif
      elseif ((j .gt. 0) .and. (k.lt.0)) then
        if (j .eq. flav) then
c--- QQbar - process bu -> bd (obtained from bu -> bd by switching
c---                               p2 and p6)
        msq(n57_6,j,k)=msq(n57_6,j,k)
     .   +sub57_6(qq)
     .   *((xn-two/xn)*m57_6(2,j,k)-m57_6(0,j,k)/xn-m57_6(1,j,k)/xn)
        msq(n67_5,j,k)=msq(n67_5,j,k)
     .   +sub67_5(qq)
     .   *((xn-two/xn)*m67_5(2,j,k)-m67_5(0,j,k)/xn-m67_5(1,j,k)/xn)
        msq(n17_5,j,k)=msq(n17_5,j,k)
     .   +sub17_5(qq)
     .   *((xn-two/xn)*m17_5(1,j,k)-m17_5(0,j,k)/xn-m17_5(2,j,k)/xn)
     .   +sub57_1(qq)
     .   *((xn-two/xn)*m17_5(1,j,k)-m17_5(0,j,k)/xn-m17_5(2,j,k)/xn)
        msq(n17_6,j,k)=msq(n17_6,j,k)
     .   +sub17_6(qq)
     .   *((xn+1d0/xn)*m17_6(0,j,k)+two*(m17_6(2,j,k)+m17_6(1,j,k))/xn)
     .   +sub67_1(qq)
     .   *((xn+1d0/xn)*m17_6(0,j,k)+two*(m17_6(2,j,k)+m17_6(1,j,k))/xn)
        msq(n27_6,j,k)=msq(n27_6,j,k)
     .   +sub27_6(qq)
     .   *((xn-two/xn)*m27_6(1,j,k)-m27_6(0,j,k)/xn-m27_6(2,j,k)/xn)
     .   +sub67_2(qq)
     .   *((xn-two/xn)*m27_6(1,j,k)-m27_6(0,j,k)/xn-m27_6(2,j,k)/xn)
        msq(n27_5,j,k)=msq(n27_5,j,k)
     .   +sub27_5(qq)
     .   *((xn+1d0/xn)*m27_5(0,j,k)+two*(m27_5(2,j,k)+m27_5(1,j,k))/xn)
     .   +sub57_2(qq)
     .   *((xn+1d0/xn)*m27_5(0,j,k)+two*(m27_5(2,j,k)+m27_5(1,j,k))/xn)
        msq(n17_2,j,k)=msq(n17_2,j,k)
     .   +sub17_2(qq)
     .   *((xn-two/xn)*m17_2(2,j,k)-m17_2(0,j,k)/xn-m17_2(1,j,k)/xn)
        msq(n27_1,j,k)=msq(n27_1,j,k)
     .   +sub27_1(qq)
     .   *((xn-two/xn)*m27_1(2,j,k)-m27_1(0,j,k)/xn-m27_1(1,j,k)/xn)
        else
c--- QQbar - process ub -> bd (obtained from ub -> bd by switching
c---                               p2 and p5)
      msq(n57_6,j,k)=msq(n57_6,j,k)
     . +sub57_6(qq)
     . *((xn-two/xn)*m57_6(1,j,k)-m57_6(0,j,k)/xn-m57_6(2,j,k)/xn)
      msq(n67_5,j,k)=msq(n67_5,j,k)
     . +sub67_5(qq)
     . *((xn-two/xn)*m67_5(1,j,k)-m67_5(0,j,k)/xn-m67_5(2,j,k)/xn)
      msq(n27_5,j,k)=msq(n27_5,j,k)
     . +sub27_5(qq)
     . *((xn-two/xn)*m27_5(2,j,k)-m27_5(0,j,k)/xn-m27_5(1,j,k)/xn)
     . +sub57_2(qq)
     . *((xn-two/xn)*m27_5(2,j,k)-m27_5(0,j,k)/xn-m27_5(1,j,k)/xn)
      msq(n27_6,j,k)=msq(n27_6,j,k)
     . +sub27_6(qq)
     . *((xn+1d0/xn)*m27_6(0,j,k)+two*(m27_6(1,j,k)+m27_6(2,j,k))/xn)
     . +sub67_2(qq)
     . *((xn+1d0/xn)*m27_6(0,j,k)+two*(m27_6(1,j,k)+m27_6(2,j,k))/xn)
      msq(n17_6,j,k)=msq(n17_6,j,k)
     . +sub17_6(qq)
     . *((xn-two/xn)*m17_6(2,j,k)-m17_6(0,j,k)/xn-m17_6(1,j,k)/xn)
     . +sub67_1(qq)
     . *((xn-two/xn)*m17_6(2,j,k)-m17_6(0,j,k)/xn-m17_6(1,j,k)/xn)
      msq(n17_5,j,k)=msq(n17_5,j,k)
     . +sub17_5(qq)
     . *((xn+1d0/xn)*m17_5(0,j,k)+two*(m17_5(1,j,k)+m17_5(2,j,k))/xn)
     . +sub57_1(qq)
     . *((xn+1d0/xn)*m17_5(0,j,k)+two*(m17_5(1,j,k)+m17_5(2,j,k))/xn)
      msq(n17_2,j,k)=msq(n17_2,j,k)
     . +sub17_2(qq)
     . *((xn-two/xn)*m17_2(1,j,k)-m17_2(0,j,k)/xn-m17_2(2,j,k)/xn)
      msq(n27_1,j,k)=msq(n27_1,j,k)
     . +sub27_1(qq)
     . *((xn-two/xn)*m27_1(1,j,k)-m27_1(0,j,k)/xn-m27_1(2,j,k)/xn)
        endif
      elseif ((j .lt. 0) .and. (k.gt.0)) then
        if (j .eq. -flav) then
c--- QbarQ - process bu -> bd (obtained from QQbar process ub -> bd
c---                               by switching p1 and p2)
      msq(n57_6,j,k)=msq(n57_6,j,k)
     . +sub57_6(qq)
     . *((xn-two/xn)*m57_6(2,j,k)-m57_6(0,j,k)/xn-m57_6(1,j,k)/xn)
      msq(n67_5,j,k)=msq(n67_5,j,k)
     . +sub67_5(qq)
     . *((xn-two/xn)*m67_5(2,j,k)-m67_5(0,j,k)/xn-m67_5(1,j,k)/xn)
      msq(n17_5,j,k)=msq(n17_5,j,k)
     . +sub17_5(qq)
     . *((xn-two/xn)*m17_5(1,j,k)-m17_5(0,j,k)/xn-m17_5(2,j,k)/xn)
     . +sub57_1(qq)
     . *((xn-two/xn)*m17_5(1,j,k)-m17_5(0,j,k)/xn-m17_5(2,j,k)/xn)
      msq(n17_6,j,k)=msq(n17_6,j,k)
     . +sub17_6(qq)
     . *((xn+1d0/xn)*m17_6(0,j,k)+two*(m17_6(2,j,k)+m17_6(1,j,k))/xn)
     . +sub67_1(qq)
     . *((xn+1d0/xn)*m17_6(0,j,k)+two*(m17_6(2,j,k)+m17_6(1,j,k))/xn)
      msq(n27_6,j,k)=msq(n27_6,j,k)
     . +sub27_6(qq)
     . *((xn-two/xn)*m27_6(1,j,k)-m27_6(0,j,k)/xn-m27_6(2,j,k)/xn)
     . +sub67_2(qq)
     . *((xn-two/xn)*m27_6(1,j,k)-m27_6(0,j,k)/xn-m27_6(2,j,k)/xn)
      msq(n27_5,j,k)=msq(n27_5,j,k)
     . +sub27_5(qq)
     . *((xn+1d0/xn)*m27_5(0,j,k)+two*(m27_5(2,j,k)+m27_5(1,j,k))/xn)
     . +sub57_2(qq)
     . *((xn+1d0/xn)*m27_5(0,j,k)+two*(m27_5(2,j,k)+m27_5(1,j,k))/xn)
      msq(n17_2,j,k)=msq(n17_2,j,k)
     . +sub17_2(qq)
     . *((xn-two/xn)*m17_2(2,j,k)-m17_2(0,j,k)/xn-m17_2(1,j,k)/xn)
      msq(n27_1,j,k)=msq(n27_1,j,k)
     . +sub27_1(qq)
     . *((xn-two/xn)*m27_1(2,j,k)-m27_1(0,j,k)/xn-m27_1(1,j,k)/xn)
        else
c--- QbarQ - process ub -> bd (obtained from QQbar process bu -> bd
c---                               by switching p1 and p2)
        msq(n57_6,j,k)=msq(n57_6,j,k)
     .   +sub57_6(qq)
     .   *((xn-two/xn)*m57_6(1,j,k)-m57_6(0,j,k)/xn-m57_6(2,j,k)/xn)
        msq(n67_5,j,k)=msq(n67_5,j,k)
     .   +sub67_5(qq)
     .   *((xn-two/xn)*m67_5(1,j,k)-m67_5(0,j,k)/xn-m67_5(2,j,k)/xn)
        msq(n27_5,j,k)=msq(n27_5,j,k)
     .   +sub27_5(qq)
     .   *((xn-two/xn)*m27_5(2,j,k)-m27_5(0,j,k)/xn-m27_5(1,j,k)/xn)
     .   +sub57_2(qq)
     .   *((xn-two/xn)*m27_5(2,j,k)-m27_5(0,j,k)/xn-m27_5(1,j,k)/xn)
        msq(n27_6,j,k)=msq(n27_6,j,k)
     .   +sub27_6(qq)
     .   *((xn+1d0/xn)*m27_6(0,j,k)+two*(m27_6(1,j,k)+m27_6(2,j,k))/xn)
     .   +sub67_2(qq)
     .   *((xn+1d0/xn)*m27_6(0,j,k)+two*(m27_6(1,j,k)+m27_6(2,j,k))/xn)
        msq(n17_6,j,k)=msq(n17_6,j,k)
     .   +sub17_6(qq)
     .   *((xn-two/xn)*m17_6(2,j,k)-m17_6(0,j,k)/xn-m17_6(1,j,k)/xn)
     .   +sub67_1(qq)
     .   *((xn-two/xn)*m17_6(2,j,k)-m17_6(0,j,k)/xn-m17_6(1,j,k)/xn)
        msq(n17_5,j,k)=msq(n17_5,j,k)
     .   +sub17_5(qq)
     .   *((xn+1d0/xn)*m17_5(0,j,k)+two*(m17_5(1,j,k)+m17_5(2,j,k))/xn)
     .   +sub57_1(qq)
     .   *((xn+1d0/xn)*m17_5(0,j,k)+two*(m17_5(1,j,k)+m17_5(2,j,k))/xn)
        msq(n17_2,j,k)=msq(n17_2,j,k)
     .   +sub17_2(qq)
     .   *((xn-two/xn)*m17_2(1,j,k)-m17_2(0,j,k)/xn-m17_2(2,j,k)/xn)
        msq(n27_1,j,k)=msq(n27_1,j,k)
     .   +sub27_1(qq)
     .   *((xn-two/xn)*m27_1(1,j,k)-m27_1(0,j,k)/xn-m27_1(2,j,k)/xn)
        endif

      elseif (j.eq.0) then
c-- G-Q
       if     (k .gt. 0) then
       if (k .eq. +flav) then
       msq(n16_2,j,k)=msq(n16_2,j,k)+sub16_2(qg)*(aveqg/aveqq)*(
     .   +(xn-one/xn)*(msq16_2(-1,k)+msq16_2(-2,k)
     .                +msq16_2(-3,k)+msq16_2(-4,k)))
       msq(n17_2,j,k)=msq(n17_2,j,k)+sub17_2(qg)*(aveqg/aveqq)*(
     .   +(xn-one/xn)*(msq17_2(+1,k)+msq17_2(+2,k)
     .                +msq17_2(+3,k)+msq17_2(+4,k)))
       else
       endif
C-- G-Qbar
        elseif (k .lt. 0) then
        if (k .eq. -flav) then
        msq(n16_2,j,k)=msq(n16_2,j,k)+sub16_2(qg)*(aveqg/aveqq)*(
     .    +(xn-one/xn)*(msq16_2(+1,k)+msq16_2(+2,k)
     .                 +msq16_2(+3,k)+msq16_2(+4,k)))
        msq(n17_2,j,k)=msq(n17_2,j,k)+sub17_2(qg)*(aveqg/aveqq)*(
     .    +(xn-one/xn)*(msq17_2(-1,k)+msq17_2(-2,k)
     .                 +msq17_2(-3,k)+msq17_2(-4,k)))
        else
        endif
        endif
      elseif (k.eq.0) then
c-- Q-G
        if     (j .gt. 0) then
        if ( j .eq. +flav) then
        msq(n26_1,j,k)=msq(n26_1,j,k)+sub26_1(qg)*(aveqg/aveqq)*(
     .    +(xn-one/xn)*(msq26_1(j,-1)+msq26_1(j,-2)
     .                 +msq26_1(j,-3)+msq26_1(j,-4)))
        msq(n27_1,j,k)=msq(n27_1,j,k)+sub27_1(qg)*(aveqg/aveqq)*(
     .    +(xn-one/xn)*(msq27_1(j,+1)+msq27_1(j,+2)
     .                 +msq27_1(j,+3)+msq27_1(j,+4)))
        else
        endif
C-- Qbar-G
        elseif (j .lt. 0) then
        if ( j .eq. -flav) then
        msq(n26_1,j,k)=msq(n26_1,j,k)+sub26_1(qg)*(aveqg/aveqq)*(
     .    +(xn-one/xn)*(msq26_1(j,+1)+msq26_1(j,+2)
     .                 +msq26_1(j,+3)+msq26_1(j,+4)))
        msq(n27_1,j,k)=msq(n27_1,j,k)+sub27_1(qg)*(aveqg/aveqq)*(
     .    +(xn-one/xn)*(msq27_1(j,-1)+msq27_1(j,-2)
     .                 +msq27_1(j,-3)+msq27_1(j,-4)))
        else
        endif
        endif
      endif

      endif
 97   continue

c--- contributions with 2 b-quarks
      if (isub .eq. 2) then
      if (j.eq.0) then
c-- G-Q
       if (kk(k) .eq. 1) then
        msq(n16_2,j,k)=msq(n16_2,j,k)+sub16_2(qg)*(aveqg/aveqq)*(
     .    +(xn-one/xn)*msq16_2(-flav,k))
        msq(n15_2,j,k)=msq(n15_2,j,k)+sub15_2(qg)*(aveqg/aveqq)*(
     .    +(xn-one/xn)*msq15_2(+flav,k))
C-- G-Qbar
        elseif (kk(k) .eq. -2) then
        msq(n16_2,j,k)=msq(n16_2,j,k)+sub16_2(qg)*(aveqg/aveqq)*(
     .    +(xn-one/xn)*msq16_2(-flav,k))
        msq(n15_2,j,k)=msq(n15_2,j,k)+sub15_2(qg)*(aveqg/aveqq)*(
     .    +(xn-one/xn)*msq15_2(+flav,k))
        endif
      elseif (k.eq.0) then
c-- Q-G
        if     (jj(j) .eq. 1) then
        msq(n26_1,j,k)=msq(n26_1,j,k)+sub26_1(qg)*(aveqg/aveqq)*(
     .    +(xn-one/xn)*msq26_1(j,+flav))
        msq(n25_1,j,k)=msq(n25_1,j,k)+sub25_1(qg)*(aveqg/aveqq)*(
     .    +(xn-one/xn)*msq25_1(j,-flav))
C-- Qbar-G
        elseif (jj(j) .eq. -2) then
        msq(n26_1,j,k)=msq(n26_1,j,k)+sub26_1(qg)*(aveqg/aveqq)*(
     .    +(xn-one/xn)*msq26_1(j,+flav))
        msq(n25_1,j,k)=msq(n25_1,j,k)+sub25_1(qg)*(aveqg/aveqq)*(
     .    +(xn-one/xn)*msq25_1(j,-flav))
        endif
      endif
      endif

      enddo
      enddo

      return
      end


text:
      SUBROUTINE PK_MISSP(KFILDO,IA,IB,NXY,MISSP,MINPK,IFIRST,
     1                    ISECOND,IMINA,NUMOCTET,SECOND,IER)
C
C        APRIL    1997   GLAHN   MODIFIED PACKXX TO HANDLE PRIMARY MISSING
C                                VALUES
C        MARCH    2000   GLAHN   REVISED FOR GRIB2;
C                                CHANGED NAME FROM PACKYY
C        JANUARY  2001   GLAHN   COMMENTS; ELIMINATED ARRAY IBXX2( );
C                                CHANGED IER = 1100 TO 920
C        NOVEMBER 2001   GLAHN   ADDED IER=0 AT BEGINNING
C        NOVEMBER 2001   GLAHN   CHANGED NUMOCTET=(NUMBIT/8)+1 TO
C                                NUMOCTET=(NUMBIT+7)/8 NEAR THE END
C        JANUARY  2002   GLAHN   ADDED TEST FOR 1ST ORDER DIFF GE
C                                2**15.; MODIFIED COMMENT FOR IER = 920
C        JANUARY  2002   GLAHN   ADDED MALLOW = 2**30+1 AND USED INSTEAD
C                                OF 999999
C        DECEMBER 2002   TAYLOR/GLAHN   CHANGED STATEMENT 225 FROM
C                                'IF(IVALUE.LT.0)NUMBIT=NUMBIT+1' TO
C                                NUMBIT=NUMBIT+1 TO ALLOW FOR SIGN
C
C        PURPOSE
C            DETERMINES WHETHER OR NOT TO USE SECOND ORDER
C            DIFFERENCES OR ORIGINAL VALUES TO PACK.  ORIGINAL VALUES
C            ARE INDICATED WHEN THE AVERAGE RANGE OF CONSECUTIVE GROUPS
C            OF SIZE MINPK OF THE SECOND ORDER DIFFERENCES IS
C            LARGER THAN THE AVERAGE RANGE OF CONSECUTIVE GROUPS OF
C            SIZE MINPK OF THE ORIGINAL VALUES.  NOTE THAT THIS
C            PROCEDURE DOES NOT IN GENERAL USE THE SAME GROUPS AS ARE
C            USED IN THE ACTUAL PACKING BECAUSE IT DOES NOT EMPLOY THE
C            LOOKBACK PROCEDURE.  THIS ALGORITHM IS RELATIVELY CHEAP
C            AND GIVES A RESULT GOOD ENOUGH FOR THE PURPOSE.  THIS
C            ROUTINE IS USED WHEN THERE CAN BE MISSING PRIMARY
C            VALUES (I.E., A MISSING VALUE INDICATOR, MISSP, CAN BE
C            IN THE DATA).
C
C        DATA SET USE
C           KFILDO - UNIT NUMBER FOR OUTPUT (PRINT) FILE. (OUTPUT)
C
C        VARIABLES
C              KFILDO = UNIT NUMBER FOR OUTPUT (PRINT) FILE.  (INPUT)
C               IA(K) = HOLDS THE NXY ORIGINAL VALUES ON INPUT
C                       (K=1,NXY).
C                       HOLDS THE NXY SECOND ORDER DIFFERENCES ON
C                       OUTPUT WHEN SECOND ORDER DIFFERENCES ARE TO BE
C                       USED.  IN THAT CASE, SECOND IS .TRUE.  SINCE
C                       THERE ARE ONLY NXY-2 VALUES, THE FIRST 2
C                       NON-MISSING
C                       VALUES ARE DUMMY.  (INPUT-OUTPUT)
C               IB(K) = WORK ARRAY (K=1,NXY).  (INTERNAL)
C                 NXY = NUMBER OF VALUES IN IA( ) ON INPUT.  ON RETURN,
C                       NXY WILL ALSO BE THE NUMBER OF VALUES IN IA( ).
C                       USED AS DIMENSION OF IA( ), IB( ), AND IWORK( ).
C                       (INPUT)
C               MISSP = PRIMARY MISSING VALUE INDICATOR.  (INPUT)
C               MINPK = INCREMENT IN WHICH RANGES WILL BE COMPUTED.
C                       (INPUT)
C              IFIRST = IFIRST IS THE FIRST ORIGINAL VALUE.  (OUTPUT)
C             ISECOND = ISECOND IS THE SECOND ORIGINAL VALUE. (OUTPUT)
C               IMINA = THE OVERALL MINIMUM VALUE OF THE SECOND ORDER
C                       SPATIAL DIFFERENCES. (OUTPUT)
C            NUMOCTET = THE MINIMUM NUMBER OF OCTETS REQUIRED TO
C                       PACK EACH OF THE VALUES REPRESENTED BY IFIRST,
C                       ISECOND, AND IMINA. (OUTPUT)
C              SECOND = TRUE (FALSE) WHEN SECOND ORDER DIFFERENCES
C                       ARE (ARE NOT) TO BE USED.  (OUTPUT)
C                 IER = CONTAINS ANY ERROR CODES GENERATED BY THIS
C                       ROUTINE
C                         0 = GOOD RETURN
C                       920 = A VALUE LARGER THAN WHAT CAN BE PACKED
C                             INTO 30 BITS HAS BEEN ENCOUNTERED.
C
C             LOCAL VARIABLES
C                AVGR = CONTAINS THE AVERAGE RANGE OF NXY VALUES
C                       IN INCREMENTS OF MINPK FOR THE ORIGINAL
C                       DATA FIELD IN IA( ).
C               AVGR2 = CONTAINS THE AVERAGE RANGE OF NXY VALUES
C                       IN INCREMENTS OF MINPK FOR THE DATA FIELD
C                       WITH SECOND ORDER DIFFERENCES IN IB( ).
C               CFEED = CONTAINS THE CHARACTER REPRESENTATION
C                       OF A PRINTER FORM FEED.
C              IAKEEP = USED IN DETERMINING SECOND ORDER
C                       DIFFERENCES AND MAINTAINING THE LOCATIONS
C                       OF MISSING VALUES.
C              ICKEEP = USED IN DETERMINING SECOND ORDER
C                       DIFFERENCES AND MAINTAINING THE LOCATIONS
C                       OF MISSING VALUES.
C               IFEED = CONTAINS THE INTEGER VALUE OF A PRINTER
C                       FORM FEED.
C              IRANGE = THE DIFFERENCE OF JMAX AND JMIN.
C             ISECPOS = THE POSITION OF THE SECOND NON-MISSING VALUE
C                       IN THE DATA FIELD.
C              IVALUE = THE VALUE OF IFIRST, ISECOND, AND IMINA
C                       THAT WILL NEED THE MOST OCTETS TO STORE.
C            IWORK(K) = WORK ARRAY. THIS ARRAY CONTAINS ALL OF THE
C                       COMPUTATIONS UNTIL WE ARE CERTAIN THAT
C                       WE WANT TO USE SECOND ORDER DIFFERENCES
C                       (K=1,NXY).  (AUTOMATIC ARRAY)
C           JMAX,JMIN = USED TO KEEP TRACK OF THE MAX AND MIN VALUES
C                       OF A PARTICULAR GROUP OF MINPK VALUES.
C              JFIRST = CONTAINS THE ABSOLUTE VALUE OF THE FIRST
C                       NON-MISSING VALUE IN THE ORIGINAL DATA
C                       FIELD.
C            JLARGEST = CONTAINS THE LARGEST OF JFIRST, JSECOND, AND
C                       JMINA FOR THE PURPOSE OF DETERMINING
C                       HOW MANY OCTETS TO STORE ALL THREE VALUES
C                       IN.
C               JMINA = CONTAINS THE ABSOLUTE VALUE OF THE MINIMUM
C                       OF THE FIELD OF SECOND ORDER DIFFERENCES.
C             JSECOND = CONTAINS THE ABSOLUTE VALUE OF THE SECOND
C                       NON-MISSING VALUE IN THE ORIGINAL DATA
C                       FIELD.
C                   K = LOOP INDEX VARIABLE.
C              KFIRST = FLAG INDICATING WHETHER THE FIRST NON-MISSING
C                       VALUE OF THE FIELD HAS BEEN FOUND.
C               KOUNT = A COUNTING/LOOPING VARIABLE.
C             KSECOND = FLAG INDICATING WHETHER THE SECOND NON-MISSING
C                       VALUE OF THE FIELD HAS BEEN FOUND.
C                SUMR = USED IN COMPUTING THE AVERAGE RANGE OF THE
C                       ORIGINAL NXY VALUES IN INCREMENTS OF MINPK.
C               LARGE = 2**15, THE LARGEST 1ST ORDER DIFFERENCE
C                       TOLERATED.  IF THE 1ST ORDER DIFFERENCE
C                       EXCEEDS THIS VALUE, IT IS LIKELY THE 2ND
C                       ORDER DIFFERENCE WILL OVERFLOW THE INTEGER
C                       COMPUTATION AND CAUSE AN UNDETECTED ERROR.
C
C        NON SYSTEM SUBROUTINES CALLED
C           NONE
C
      PARAMETER (MALLOW=2**30+1)
      PARAMETER (LARGE=2**15)
C
      CHARACTER*1 CFEED
      LOGICAL SECOND,KFIRST,KSECOND
C
      DIMENSION IA(NXY),IB(NXY),IWORK(NXY)
C
      DATA IFEED/12/
C
      IER=0
C
C        INITIALIZE THE FORM FEED CHARACTER.
      CFEED=CHAR(IFEED)
C
C        INITIALIZE IFIRST, ISECOND, IB(1), AND ICKEEP.
      IFIRST=IA(1)
      ISECOND=MISSP
      IB(1)=0
      IF(IA(1).EQ.MISSP)IB(1)=MISSP
      ICKEEP=IA(1)
C
C        COMPUTE FIRST ORDER DIFFERENCES.
C
      DO 120 K=2,NXY
C
         IF(ICKEEP.EQ.MISSP)THEN
            IB(K)=MISSP
            ICKEEP=IA(K)
         ELSE
C
            IF(IA(K).EQ.MISSP)THEN
               IB(K)=MISSP
            ELSE
               IB(K)=IA(K)-ICKEEP
C
               IF(IB(K).GT.LARGE)THEN
                  SECOND=.FALSE.
C                 WRITE(KFILDO,119)
C119              FORMAT(/' 1ST ORDER DIFFERENCE EXCEEDS 2**15.',
C    1                    '  SECOND ORDER DIFFERENCES ARE NOT PAKCED.')
                  GO TO 900
               ENDIF
C
               ICKEEP=IA(K)
            ENDIF
C
         ENDIF
C
         IF(IFIRST.EQ.MISSP)THEN
C
            IF(IA(K).NE.MISSP)THEN
               IFIRST=IA(K)
               IB(K)=0
            ENDIF
C
         ELSE IF(ISECOND.EQ.MISSP)THEN
C
            IF(IA(K).NE.MISSP)THEN
               ISECOND=IA(K)
            ENDIF
C
         ENDIF
C
 120  CONTINUE
C
C***D     WRITE(KFILDO,121)(IA(K),K=1,10000)
C***D     WRITE(KFILDO,121)(IA(K),K=268351,268500)
C***D121  FORMAT(/' ORIGINAL SCALED VALUES'/(' '25I4))
C
C***D     WRITE(KFILDO,1210)IFIRST,IFOD
C***D1210 FORMAT(/' IFIRST ='I11,'     IFOD ='I11)
C***D     WRITE(KFILDO,122)(IB(K),K=268351,268500)
C***D122  FORMAT(/ 'FIRST ORDER DIFFERENCES'/(' '25I4))
C
C        COMPUTE SECOND ORDER DIFFERENCES.
C
      IAKEEP=MISSP
      KFIRST=.TRUE.
      KSECOND=.TRUE.
C
      DO 130 K=1,NXY
C
         IF(KFIRST)THEN
            IWORK(K)=IB(K)
            IF(IB(K).NE.MISSP)KFIRST=.FALSE.
         ELSE
C
            IF(IAKEEP.EQ.MISSP)THEN
C
               IF(IB(K).NE.MISSP)THEN
                  IAKEEP=IB(K)
                  IF(KSECOND)THEN
                     IB(K)=0
                     KSECOND=.FALSE.
                     ISECPOS=K
                  ENDIF
                  IWORK(K)=IB(K)
               ELSE
                  IWORK(K)=MISSP
               ENDIF
C
            ELSE
C
               IF(IB(K).EQ.MISSP)THEN
                  IWORK(K)=MISSP
               ELSE
                  IWORK(K)=IB(K)-IAKEEP
                  IMINA=IWORK(K)
                  IAKEEP=IB(K)
               ENDIF
C
            ENDIF
C
         ENDIF
C
 130  CONTINUE
C
C***D     WRITE(KFILDO,136)CFEED
C***D136  FORMAT(1A/' *********************'
C***D    1          /' 2ND ORDER DIFFERENCES'
C***D    2          /' *********************')
C***D     WRITE(KFILDO,137) (IWORK(K),K=268351,268500)
C***D137  FORMAT(/(' '25I4))
C
C        COMPUTE THE MINIMUM OF THE FIELD OF SECOND ORDER
C        DIFFERENCES AND SUBTRACT IT FROM OF EACH OF THE
C        VALUES IN THE FIELD.
C
      DO 132 K=ISECPOS+1,NXY
         IF(IWORK(K).EQ.MISSP)CYCLE
         IF(IWORK(K).LT.IMINA)IMINA=IWORK(K)
 132  CONTINUE
C
      DO 134 K=ISECPOS+1,NXY
         IF(IWORK(K).EQ.MISSP)CYCLE
         IWORK(K)=IWORK(K)-IMINA
 134  CONTINUE
C
C***D     WRITE(KFILDO,131)(IWORK(K),K=268351,268500)
C***D131  FORMAT(/' SECOND ORDER DIFFERENCES'/(' '25I4))
C
C        COMPUTE AVERAGE RANGE OF NXY ORIGINAL VALUES IN INCREMENTS OF
C        MINPK.
C
      SUMR=0
      KOUNT=0
C
      DO 140 K=1,NXY,MINPK
         JMIN=MALLOW
         JMAX=-MALLOW
         IF(K+MINPK-1.GT.NXY)GO TO 140
C           THE LAST GROUP MAY BE VERY SMALL AND NOT BE REPRESENTATIVE
C           OF THE RANGE.
C
         DO 135 J=K,K+MINPK-1
            IF(IA(J).EQ.MISSP)CYCLE
            IF(IA(J).GT.JMAX)JMAX=IA(J)
            IF(IA(J).LT.JMIN)JMIN=IA(J)
 135     CONTINUE
C
         KOUNT=KOUNT+1
C
         IF(JMIN.EQ.MALLOW)THEN
            IRANGE=0
         ELSE
            IRANGE=JMAX-JMIN
            SUMR=SUMR+IRANGE
         ENDIF
C
 140  CONTINUE
C
      AVGR=99999.
      IF(KOUNT.NE.0)AVGR=SUMR/KOUNT
C
C        COMPUTE AVERAGE RANGE OF NXY-4 2ND ORDER VALUES IN
C        INCREMENTS OF MINPK.   DON'T USE THE FIRST 2 VALUES,
C        BECAUSE THEY COULD BE BASED ON MISSING VALUES.
C
      SUMR=0
      KOUNT=0
C
      DO 150 K=3,NXY-2,MINPK
         JMIN=MALLOW
         JMAX=-MALLOW
         IF(K+MINPK-1.GT.NXY-2)GO TO 150
C        THE LAST GROUP MAY BE VERY SMALL AND NOT BE REPRESENTATIVE OF
C        THE RANGE.
C
         DO 145 J=K,K+MINPK-1
            IF(IB(J).EQ.MISSP)GO TO 145
            IF(IB(J).GT.JMAX)JMAX=IB(J)
            IF(IB(J).LT.JMIN)JMIN=IB(J)
 145     CONTINUE
C
         KOUNT=KOUNT+1
C
         IF(JMIN.EQ.MALLOW)THEN
            IRANGE=0
         ELSE
            IRANGE=JMAX-JMIN
            SUMR=SUMR+IRANGE
         ENDIF
C
 150  CONTINUE
C
      AVGR2=99999.
      IF(KOUNT.NE.0)AVGR2=SUMR/KOUNT
C
      SECOND=.FALSE.
C     WRITE(KFILDO,155)AVGR,AVGR2
C155  FORMAT(/' AVERAGE RANGE OF ORIGINAL SCALED VALUES =  'F10.2/
C    1        ' AVERAGE RANGE OF SECOND ORDER DIFFERENCES ='F10.2)
      IF(AVGR2.GE.AVGR)GO TO 900
C        SECOND ORDER DIFFERENCES NOT TO BE USED.
C
C        TRANSFER IWORK( ) TO IA( ).
      KOUNT=0
C
      DO 180 K=1,NXY
         IA(K)=IWORK(K)
 180  CONTINUE
C
C***D     WRITE(KFILDO,185)(IA(K),K=1,100)
C***D185  FORMAT(/' SECOND ORDER DIFFERENCES TO BE PACKED'/(' '10I11))
C
C           SINCE WE ARE DOING SECOND ORDER DIFFERENCES
C           TAKE THE TIME HERE TO DETERMINE THE MINIMUM
C           NUMBER OF BYTES (NOT BITS) NEED TO CONTAIN
C           THE FIRST ORIGINAL VALUE, THE SECOND ORIGINAL
C           VALUE, AND THE MINIMUM OF THE FIELD OF SECOND
C           ORDER DIFFERENCES.
C
      NUMBIT=0
      JMINA=ABS(IMINA)
      JFIRST=ABS(IFIRST)
      JSECOND=ABS(ISECOND)
C
C           FIND THE LARGEST ABSOLUTE VALUE OF THE
C           THE THREE VARIABLES IMINA, IFIRST, AND
C           ISECOND.
C
      IF(JMINA.GT.JFIRST)THEN
         JLARGEST=JMINA
         IVALUE=IMINA
      ELSE
         JLARGEST=JFIRST
         IVALUE=IFIRST
      ENDIF
C
      IF(JLARGEST.LT.JSECOND)THEN
         JLARGEST=JSECOND
         IVALUE=ISECOND
      ENDIF
C
      IBXX2=2
C
      DO 220 K=1,30
         NUMBIT=K
         IF(JLARGEST.LT.IBXX2)GO TO 225
         IBXX2=IBXX2*2
 220  CONTINUE
C
C        DROP THROUGH HERE MEANS THE NUMBER OF BITS
C        EXCEEDS 30 AND THE VALUE TO PACK EXCEEDS 2**30.
C        2**31 IS TOO LARGE A VALUE TO PACK.  2**30
C        IS USED AS THE LIMIT THROUGHOUT THE PACKER,
C        ALTHOUGH A VALUE OF 2**31-1 COULD PROBABLY
C        BE ACCOMMODATED.
      IER=920
C        THIS IS A FATAL ERROR.
      GO TO 900
C
 225  NUMBIT=NUMBIT+1
C        ALLOW FOR SIGN BIT.
C
C        FIND THE SMALLEST NUMBER OF OCTETS
C        THAT WILL CONTAIN THE VALUES.
      NUMOCTET=(NUMBIT+7)/8
      SECOND=.TRUE.
C     WRITE(KFILDO,899)JMINA,JFIRST,JSECOND,JLARGEST,
C    1                 IVALUE,NUMBIT,NUMOCTET
C899  FORMAT(/' AT 899 IN PK_MISSP--JMINA,JFIRST,JSECOND,JLARGEST,',
C    1        'IVALUE,NUMBIT,NUMOCTET'/20X,7I10)
 900  RETURN
      END

ext: f

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
program gen_be_ensrf
!
!----------------------------------------------------------------------
!  Purpose : Perform an Ensemble Square Root Filter (EnSRF) assimilation
!  using WRF forecast data.
!
!  Owner: Dale Barker (NCAR/MMM)
!  Please acknowledge author/institute in work that uses this code.
!
!----------------------------------------------------------------------

#ifdef crayx1
#define iargc ipxfargc
#endif

   use da_control, only : stdout,stderr,filename_len
   use da_gen_be, only : da_stage0_initialize, da_get_field, da_get_trh

   implicit none

   type info_type
       character*40   :: id
       character*40   :: platform
       integer        :: plat
       integer        :: variable
   end type info_type

   type loc_type
       real           :: x
       real           :: y
       real           :: z
   end type loc_type

   type ob_type
      type (info_type):: info
      type (loc_type) :: loc
      real            :: yo
      real            :: omb_mean
      real            :: oma_mean
      real            :: omb_rms
      real            :: oma_rms
      real            :: sigma_o
      real            :: sigma_b
   end type ob_type

   integer, parameter    :: nv3d = 4                  ! #3D variables (u, v, T, q).
   integer, parameter    :: nv2d = 1                  ! #2D variables (ps).
   integer, parameter    :: unit = 100                ! Unit number.

   character (len=filename_len)   :: filestub                  ! General filename stub.
   character (len=filename_len)   :: input_file                ! Input file.
   character (len=filename_len)   :: output_file               ! Output file.
   character (len=10)    :: var                       ! Variable to search for.
   character (len=3)     :: ce                        ! Member index -> character.

   integer               :: num_members               ! Ensemble size.
   integer               :: num_obs                   ! Number of observations.
   integer               :: ni                        ! 1st dimension size.
   integer               :: niu                       ! 1st dimension size (u)
   integer               :: nj                        ! 2nd dimension size.
   integer               :: njv                       ! 2nd dimension size (v)
   integer               :: nk                        ! 3rd dimension size.
   integer               :: nv                        ! Total number of 2D/3D variables.
   integer               :: nij, nijk, nijkv          ! Dimensions.
   integer               :: o, member, i, j, k, v     ! Loop counters.
   integer               :: imin, imax, jmin, jmax    ! Min/max i,j for localization.
   integer               :: ijkv                      !
   integer               :: index                     !
   integer               :: xend
   real                  :: num_members_inv
   real                  :: num_members_inv1
   real                  :: cov_inf_fac               ! Covariance inflation factor.
   real                  :: cov_loc_rad_m             ! Covariance localization radius (m).
   real                  :: cov_loc_rad               ! Covariance localization radius (gridpts).
   real                  :: ds                        ! Grid resolution (m).
   real                  :: x1, y1                    ! Grid location of ob.
   real                  :: sigma_o2, sigma_b2
   real                  :: sigma_hh_inv
   real                  :: sigma_xh
   real                  :: beta
   real                  :: o_minus_b                 ! Observation minus forecast.
   real                  :: kalman_gain
   real                  :: x_mean
   real                  :: h_mean
   real                  :: hb_mean
   real                  :: ha_mean
   real                  :: sum_f, sum_a
   real                  :: hf, ha

   type (ob_type),pointer:: ob(:)

   character(len=10),allocatable :: varr(:)          ! Variable name.
   integer, allocatable  :: nkk(:)              ! Vertical dimension of field.
   integer, allocatable  :: xstart(:)           ! Starting position of variable..

   real, allocatable     :: x(:,:), xf(:,:)
   ! real, allocatable     :: xf_cif(:,:)
   real, allocatable     :: h(:), h_devn(:)
   real, allocatable     :: x_devn(:)

   real, allocatable     :: uc(:,:)                   ! u-wind (C grid).
   real, allocatable     :: vc(:,:)                   ! v-wind (C grid).
   real, allocatable     :: field(:,:)                ! Input field.
   real, allocatable     :: dummy(:,:)                ! Dummy.
   real, allocatable     :: corr(:,:)                 ! Correlation.
   real, allocatable     :: xf_mean(:)                ! Prior mean.
   real, allocatable     :: xa_mean(:)                ! Posterior mean.

   namelist / gen_be_ensrf_nl / filestub, num_members, &
                                cov_inf_fac, cov_loc_rad_m

   stderr = 0
   stdout = 6

!---------------------------------------------------------------------------------------------
   write(6,'(/a)')' [1] Initialize information.'
!---------------------------------------------------------------------------------------------

   filestub = 'test'
   num_members = 56
   cov_inf_fac = 1.0
   cov_loc_rad_m = 1500000.0

   open(unit=unit, file='gen_be_ensrf_nl.nl', &
        form='formatted', status='old', action='read')
   read(unit, gen_be_ensrf_nl)
   close(unit)

   write(6,'(a,a)')'   Filestub = ', trim(filestub)
   write(6,'(a,i4)')'   Number of ensemble members = ', num_members
   write(6,'(a,f16.8)')'   Covariance Inflation Factor = ', cov_inf_fac
   write(6,'(a,f16.8)')'   Covariance Localization Radius (m) = ', cov_loc_rad_m

   num_members_inv = 1.0 / real(num_members)
   num_members_inv1 = 1.0 / real(num_members-1)

!---------------------------------------------------------------------------------------------
   write(6,'(/a)')' [2] Set up data dimensions and allocate arrays:'
!---------------------------------------------------------------------------------------------

!  Get grid dimensions from first T field:
   input_file = trim(filestub)//'.e001'
   var = 'T'
   call da_stage0_initialize( input_file, var, ni, nj, nk, ds )
   niu = ni+1 ! u i dimension is 1 larger.
   njv = nj+1 ! v j dimension is 1 larger.
   cov_loc_rad = cov_loc_rad_m / ds
   nij = ni * nj
   nijk = nij * nk
   nijkv = nv3d * nijk + nv2d * nij ! #3D + #2D variables.
   nv = nv3d + nv2d
print *, ni, nj, nk, nv, nij, nijk, nv3d, nv2d, nijkv
!  Allocate arrays:
   allocate ( varr(1:nv) )
   allocate ( nkk(1:nv) )
   allocate ( xstart(1:nv) )
   allocate( x(1:nijkv,1:num_members) )
   allocate( xf(1:nijkv,1:num_members) )
   ! allocate( xf_cif(1:nijkv,1:num_members) )
   allocate( h(1:num_members) )
   allocate( h_devn(1:num_members) )
   allocate( x_devn(1:num_members) )
   do v = 1, nv3d
      nkk(v) = nk
   end do
   do v = nv3d + 1, nv
      nkk(v) = 1
   end do
   allocate( field(1:ni,1:nj) )
   allocate( uc(1:niu,1:nj) )
   allocate( vc(1:ni,1:njv) )
   allocate( dummy(1:ni,1:nj) )
   allocate( corr(1:ni,1:nj) )

!  Hardwired:
   varr(1) = 'U'
   varr(2) = 'V'
   varr(3) = 'T'
   varr(4) = 'QVAPOR'
   varr(5) = 'PSFC'

!---------------------------------------------------------------------------------------------
   write(6,'(/a)')' [3] Read observations.'
!---------------------------------------------------------------------------------------------

   input_file = 'observations'
   open(unit, file = input_file, status='old')
!write(56,*)1
!      write(56,'(2a10,i6,5f16.8)')'test', 'dale', 3, &
!            20.0, 20.0, 12.0, 250.0, 1.0
!stop
   read(unit,*)num_obs
   write(6,'(a,i10)')'   Number of observations = ', num_obs
   allocate( ob(1:num_obs) )

   do o = 1, num_obs
      read(unit,'(2a10,i6,5f16.8)')ob(o) % info % id, ob(o) % info % platform, &
                ob(o) % info % variable, &
                ob(o) % loc % x, ob(o) % loc % y ,ob(o) % loc % z, &
                ob(o) % yo, ob(o) % sigma_o
   end do

!---------------------------------------------------------------------------------------------
   write(6,'(/a)')' [4] Extract necessary fields from WRF ensemble forecasts.'
!---------------------------------------------------------------------------------------------

   xstart(1) = 1
   do v = 2, nv
      xstart(v) = xstart(v-1) + ( ni * nj * nkk(v-1) )
   end do

   do member = 1, num_members

      write(UNIT=ce,FMT='(i3.3)')member
      input_file = trim(filestub)//'.e'//ce

      do v = 1, nv
         do k = 1, nkk(v)

            if ( varr(v) == 'U' ) then
               call da_get_field( input_file, varr(v), 3, niu, nj, nk, k, uc )
               do j = 1, nj
                  do i = 1, ni
                     field(i,j) = 0.5 * ( uc(i,j) + uc(i+1,j) )
                  end do
               end do
            end if

            if ( varr(v) == 'V' ) then
               call da_get_field( input_file, varr(v), 3, ni, njv, nk, k, vc )
               do j = 1, nj
                  do i = 1, ni
                     field(i,j) = 0.5 * ( vc(i,j) + vc(i,j+1) )
                  end do
               end do
            end if

            if ( varr(v) == "T" ) then
!              Read theta, and convert to temperature:
               call da_get_trh( input_file, ni, nj, nk, k, field, dummy )
            end if

!           Read mixing ratio, and convert to specific humidity:
            if ( varr(v) == 'QVAPOR' ) then
               call da_get_field( input_file, varr(v), 3, ni, nj, nk, k, field )
               field(:,:) = field(:,:) / ( 1.0 + field(:,:) )
            end if

!           Read surface pressure:
            if ( varr(v) == 'PSFC' ) then
               call da_get_field( input_file, varr(v), 2, ni, nj, nk, k, field )
            end if

!           Fill 4D array:
            index = xstart(v) + (k-1) * nij
            do j = 1, nj
               do i = 1, ni
                  xf(index,member) = field(i,j)
                  index = index + 1
               end do
            end do
         end do ! k
      end do ! v
   end do !member

!  Initialize analysis as first guess:
   x(1:nijkv,1:num_members) = xf(1:nijkv,1:num_members)

!  Perform initial covariance inflation (to boost input be):
   call da_cov_inflation( nijkv, num_members, cov_inf_fac, x )
!   xf_cif(1:nijkv,1:num_members) = x(1:nijkv,1:num_members) ! Diagnostic

   deallocate( uc )
   deallocate( vc )
   deallocate( field )
   deallocate( dummy )

!---------------------------------------------------------------------------------------------
   write(6,'(/a)')' [5] Perform EnSRF:'
!---------------------------------------------------------------------------------------------

   do o = 1, num_obs
      write(6,'(2(a,i8))')'     Assimilating observation ', o, ' of ', num_obs
      x1 = ob(o) % loc % x
      y1 = ob(o) % loc % y
      sigma_o2 = ob(o) % sigma_o * ob(o) % sigma_o

!     Perform observation operator:

      do member = 1, num_members
         call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), x(:,member), h(member) )
      end do

!     Mean and deviations of H(x):
      h_mean = sum(h(1:num_members)) * num_members_inv
      h_devn(1:num_members) = h(1:num_members) - h_mean

!     Covariance H Pf H^T:
      sigma_b2 = sum(h_devn(1:num_members) * h_devn(1:num_members)) * num_members_inv1

      ob(o) % sigma_b = sqrt(sigma_b2)

!     Store 1/(H Pf HT + R):
      sigma_hh_inv = 1.0 / ( sigma_b2 + sigma_o2 )

!     Calculate EnSRF deviation Kalman Gain scaling:
      beta = 1.0 / ( 1.0 + sqrt( sigma_o2 * sigma_hh_inv ) )

!     Ensemble mean innovation vector:
!!!DALE Uncomment for pseudo observation
!      print *, 'Running pseudo observation O-B=1.0:'
!      ob(o) % yo = h_mean + 1.0
!!!DALE
      o_minus_b = ob(o) % yo - h_mean
      write(6,'(a,2f15.5)')' Observation increment = ', o_minus_b
      write(6,'(a,2f15.5)')' Observation/forecast error variance = ', sigma_o2, sigma_b2

!-----------------------------------------------------------------------------
!     Calculate local support:
!-----------------------------------------------------------------------------

      imin = max(nint(ob(o) % loc % x - cov_loc_rad),1)
      imax = min(nint(ob(o) % loc % x + cov_loc_rad),ni)
      jmin = max(nint(ob(o) % loc % y - cov_loc_rad),1)
      jmax = min(nint(ob(o) % loc % y + cov_loc_rad),nj)

!      write(6,'(a,2f8.2)')'   Ob location(x/y) = ', ob(o) % loc % x, ob(o) % loc % y
!      write(6,'(a,4i5)')'   Min/max i/j = ', imin, imax, jmin, jmax

      do j = jmin, jmax
         do i = imin, imax
            call compact_support( x1, real(i), y1, real(j), cov_loc_rad, corr(i,j) )
         end do
      end do

      do v = 1, nv
         do k = 1, nkk(v)
            do j = jmin, jmax
               do i = imin, imax
                  ijkv = (v-1) * nijk + (k-1) * nij + (j-1) * ni + i

                  x_mean = sum(x(ijkv,1:num_members)) * num_members_inv

                  x_devn(1:num_members) = x(ijkv,1:num_members) - x_mean

!                 Calculate PfHT:
                  sigma_xh = sum( x_devn(1:num_members) * h_devn(1:num_members) ) * &
                             num_members_inv1

!                 Apply covariance localization:
                  sigma_xh = sigma_xh * corr(i,j)

                  kalman_gain = sigma_xh * sigma_hh_inv

                  x_mean = x_mean + kalman_gain * o_minus_b

                  x_devn(1:num_members) = x_devn(1:num_members) - &
                                          beta * kalman_gain * h_devn(1:num_members)

                  x(ijkv,1:num_members) = x_mean + x_devn(1:num_members)

               end do
            end do
         end do
      end do
   end do

!---------------------------------------------------------------------------------------------
   write(6,'(/a)')' [4] Calculate diagnostics:'
!---------------------------------------------------------------------------------------------

!  Calculate mean analysis over all grid-points (diagnostic only):

   allocate( xf_mean(1:nijkv) )
   allocate( xa_mean(1:nijkv) )

   do ijkv = 1, nijkv
      xf_mean(ijkv) = sum(xf(ijkv,1:num_members)) * num_members_inv
      xa_mean(ijkv) = sum(x(ijkv,1:num_members)) * num_members_inv
   end do

!  Calculate grid-point diagnostics:

!   call da_get_grid_stats( num_members, ni, nj, nk, nv, nijkv, nkk, varr, xf_cif, x )
   call da_get_grid_stats( num_members, ni, nj, nk, nv, nijkv, nkk, varr, xf, x )

!  Calculate observation diagnostics:

   write(53,*)num_obs
   do o = 1, num_obs

!     Ensemble mean-observation diagnostics:
      call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), xf_mean, hb_mean )
      call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), xa_mean, ha_mean )

      ob(o) % omb_mean = ob(o) % yo - hb_mean
      ob(o) % oma_mean = ob(o) % yo - ha_mean

!     Ensemble spread-observation diagnostics:
      sum_f = 0.0
      sum_a = 0.0
      do member = 1, num_members
!         call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), xf_cif(:,member), hf )
         call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), xf(:,member), hf )
         call da_obs_operator( nv, ni, nj, nkk, xstart, ob(o), x(:,member), ha )
         sum_f = sum_f + ( ob(o) % yo - hf )**2
         sum_a = sum_a + ( ob(o) % yo - ha )**2
      end do

      ob(o) % omb_rms = sqrt( sum_f * num_members_inv )
      ob(o) % oma_rms = sqrt( sum_a * num_members_inv )

      write(53,'(2a10,i8,10f16.8)')ob(o) % info % id, ob(o) % info % platform, &
                ob(o) % info % variable, &
                ob(o) % loc % x, ob(o) % loc % y ,ob(o) % loc % z, &
                ob(o) % yo, ob(o) % omb_mean, ob(o) % oma_mean, &
                ob(o) % omb_rms, ob(o) % oma_rms, &
                ob(o) % sigma_o, ob(o) % sigma_b
   end do

!---------------------------------------------------------------------------------------------
   write(6,'(/a)')' [4] Output EnSRF analysis ensemble:'
!---------------------------------------------------------------------------------------------

   do member = 1, num_members
      write(6,'(a,i4)')'   Writing ensemble member increments for member', member
      write(UNIT=ce,FMT='(i3.3)')member

      do v = 1, nv

!        Output prior ensemble forecasts:
         output_file = trim(varr(v))//'/'//trim(varr(v))//'.prior.e'//trim(ce)
         open(unit, file = output_file, form='unformatted')
         write(unit)ni, nj, nkk(v)
         xend = xstart(v) + nij * nkk(v) - 1
         write(unit)xf(xstart(v):xend,member)
         close(unit)

!        Output posterior ensemble forecasts:
         output_file = trim(varr(v))//'/'//trim(varr(v))//'.posterior.e'//trim(ce)
         open(unit, file = output_file, form='unformatted')
         write(unit)ni, nj, nkk(v)
         xend = xstart(v) + nij * nkk(v) - 1
         write(unit)x(xstart(v):xend,member)
         close(unit)
      end do
   end do

   do v = 1, nv

!     Output prior ensemble forecasts:
      output_file = trim(varr(v))//'/'//trim(varr(v))//'.prior.mean'
      open(unit, file = output_file, form='unformatted')
      write(unit)ni, nj, nkk(v)
      xend = xstart(v) + nij * nkk(v) - 1
      write(unit)xf_mean(xstart(v):xend)
      close(unit)

!     Output posterior ensemble forecasts:
      output_file = trim(varr(v))//'/'//trim(varr(v))//'.posterior.mean'
      open(unit, file = output_file, form='unformatted')
      write(unit)ni, nj, nkk(v)
      xend = xstart(v) + nij * nkk(v) - 1
      write(unit)xa_mean(xstart(v):xend)
      close(unit)
   end do

   deallocate( ob )

contains

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

subroutine da_obs_operator( nv, ni, nj, nkk, xstart, ob, x, h )

   implicit none

   integer, intent(in)              :: nv, ni, nj
   integer, intent(in)              :: nkk(1:nv)
   integer, intent(in)              :: xstart(1:nv)
   type (ob_type), intent(in)       :: ob
   ! WHY? was this a target?
   ! real, intent(in), target         :: x(:)
   real, intent(in)                 :: x(:)
   real, intent(out)                :: h

   integer                          :: v, nk
   integer                          :: start
   real                             :: ob_x, ob_y, ob_z ! Observation location.

   v = ob % info % variable
   nk = nkk(v)
   start = xstart( v )

!  Perform 3D interpolation:

   ob_x = ob % loc % x
   ob_y = ob % loc % y
   ob_z = ob % loc % z

   call da_interpolate_3d( ni, nj, nk, ob % loc % x, ob % loc % y, ob % loc % z, &
                           x(start:start + ni * nj * nk ), h )

end subroutine da_obs_operator

subroutine da_interpolate_3d( ni, nj, nk, x, y, z, field, h )

   implicit none
   integer, intent(in)   :: ni, nj, nk
   real, intent(in)      :: x, y, z      !  Grid location of point.
   real, intent(in)      :: field(1:ni,1:nj,1:nk) ! Field to interpolate.
   real, intent(out)     :: h            ! Interpolated value.

   integer               :: ii, jj, kk   ! Grid locators.
   real                  :: dx, dy, dz   ! Interpolation weights.
   real                  :: h1, h2       ! Interpolation values.
   real                  :: h_below      ! 2D interp values.
   real                  :: h_above      ! 2D interp values.

   ii = int(x)
   jj = int(y)
   kk = int(z)

   dx = x - real(ii)
   dy = y - real(jj)
   dz = z - real(kk)

!  2D interpolation on layer below:
   h1 = ( 1.0 - dx ) * field(ii,jj,kk)   + dx * field(ii+1,jj,kk)
   h2 = ( 1.0 - dx ) * field(ii,jj+1,kk) + dx * field(ii+1,jj+1,kk)
   h_below = ( 1.0 - dy ) * h1 + dy * h2

!  2D interpolation on layer above:
   h1 = ( 1.0 - dx ) * field(ii,jj,kk+1)   + dx * field(ii+1,jj,kk+1)
   h2 = ( 1.0 - dx ) * field(ii,jj+1,kk+1) + dx * field(ii+1,jj+1,kk+1)
   h_above = ( 1.0 - dy ) * h1 + dy * h2

!  Interpolation in vertical:
   h = ( 1.0 - dz ) * h_below + dz * h_above

end subroutine da_interpolate_3d

subroutine compact_support( x1, x2, y1, y2, cov_loc_rad, corr )

!  Compact support according to (4.10) of Gaspari and Cohn (1999).

!  Assumes r >=0.

   implicit none

   real, intent(in)     :: x1, x2         ! x values of two points.
   real, intent(in)     :: y1, y2         ! y values of two points.
   real, intent(in)     :: cov_loc_rad    ! Cut-off correlation > 2c.
   real, intent(out)    :: corr           ! Compactly-supported correlation.

   real                 :: z              ! Distance between points.
   real                 :: r              ! Ratio.
   real                 :: r2, r3, r4, r5 ! Powers of r.

   z = sqrt( ( x2 - x1 )**2 + ( y2 - y1)**2 )
   r = z / cov_loc_rad

   r2 = r * r
   r3 = r2 * r
   r4 = r3 * r
   r5 = r4 * r

   if ( r <= 1.0 ) then
      corr = -0.25 * r5 + 0.5 * r4 + 0.625 * r3 - 5.0 * r2 / 3.0 + 1
   else if ( r > 1.0 .and. r < 2.0 ) then
      corr = r5 / 12.0 - 0.5 * r4 + 0.625 * r3 + 5.0 * r2 / 3.0 - &
             5.0 * r + 4.0 - 2.0 / ( 3.0 * r )

   else if ( r >= 2.0 ) then
      corr = 0.0
   end if

end subroutine compact_support

subroutine da_get_grid_stats( num_members, ni, nj, nk, nv, nijkv, nkk, varr, xf, xa )

   implicit none

   integer, intent(in)   :: num_members
   integer, intent(in)   :: ni
   integer, intent(in)   :: nj
   integer, intent(in)   :: nk
   integer, intent(in)   :: nv
   integer, intent(in)   :: nijkv
   integer, intent(in)   :: nkk(1:nv)
   character*10,intent(in):: varr(1:nv)
   real, intent(in)      :: xf(1:nijkv,1:num_members)
   real, intent(in)      :: xa(1:nijkv,1:num_members)

   character*10          :: name
   integer               :: v, k, j, i, ijkv
   real                  :: num_members_inv, nij_inv
   real                  :: mnsq_f, mnsq_a
   real                  :: domain_mean_f, domain_mean_a
   real                  :: domain_stdv_f, domain_stdv_a
   real                  :: mean_f(1:ni,1:nj), mean_a(1:ni,1:nj)
   real                  :: stdv_f(1:ni,1:nj), stdv_a(1:ni,1:nj)

   num_members_inv = 1.0 / num_members
   nij_inv = 1.0 / real ( ni * nj )

   write(6,'(a)')' Variable, Level,   Mean_F,   Mean_A,   StDv_F,   StdV_A'

   ijkv = 0
   do v = 1, nv
      name = varr(v)
      do k = 1, nkk(v)
         do j = 1, nj
            do i = 1, ni
               ijkv = ijkv + 1
               mean_f(i,j) = sum(xf(ijkv,1:num_members)) * num_members_inv
               mean_a(i,j) = sum(xa(ijkv,1:num_members)) * num_members_inv
               mnsq_f      = sum(xf(ijkv,1:num_members)**2) * num_members_inv
               mnsq_a      = sum(xa(ijkv,1:num_members)**2) * num_members_inv
               stdv_f(i,j) = mnsq_f - mean_f(i,j) * mean_f(i,j)
               stdv_a(i,j) = mnsq_a - mean_a(i,j) * mean_a(i,j)
               stdv_f(i,j) = 0.0
               if ( stdv_f(i,j) > 0.0 ) stdv_f(i,j) = sqrt(stdv_f(i,j))
               stdv_a(i,j) = 0.0
               if ( stdv_a(i,j) > 0.0 ) stdv_a(i,j) = sqrt(stdv_a(i,j))
            end do
         end do

         domain_mean_f = sum(mean_f(:,:)) * nij_inv
         domain_mean_a = sum(mean_a(:,:)) * nij_inv
         domain_stdv_f = sum(stdv_f(:,:)) * nij_inv
         domain_stdv_a = sum(stdv_a(:,:)) * nij_inv

         write(6,'(a,i6,4f16.8)')trim(name), k, domain_mean_f, domain_mean_a, &
                                 domain_stdv_f, domain_stdv_a
      end do
   end do

end subroutine da_get_grid_stats

subroutine da_cov_inflation( nijkv, num_members, cov_inf_fac, x )

   implicit none

   integer, intent(in) :: nijkv
   integer, intent(in) :: num_members
   real, intent(in)    :: cov_inf_fac
   real, intent(inout) :: x(1:nijkv,1:num_members)

   integer             :: ijkv
   real                :: num_members_inv
   real                :: x_devn(1:num_members)

   num_members_inv = 1.0 / real(num_members)

   do ijkv = 1, nijkv
      x_mean = sum(x(ijkv,1:num_members)) * num_members_inv
      x_devn(1:num_members) = cov_inf_fac * ( x(ijkv,1:num_members) - x_mean )
      x(ijkv,1:num_members) = x_mean + x_devn(1:num_members)
   end do

end subroutine da_cov_inflation

end program gen_be_ensrf

max_stars_repo_name: matzegoebel/WRF-fluxavg

max_issues_repo_name: teb-model/wrf-teb

max_forks_repo_name: teb-model/wrf-teb


text:
subroutine colinearity(axis, origin, coords, n, x, nCon)
  implicit none
  real, dimension(3)::origin, axis
  real, dimension(n,3)::coords
  real, dimension(nCon)::x
  integer:: n, nCon

  real, dimension(n,3)::dirVec
  real, dimension(3)::resultDir
  integer:: i,j

  !Compute the direction from each point to the origin
  do i = 1,n
     do j = 1,3
        dirVec(i,j) = origin(j)-coords(i,j)
     end do
  end do

  ! compute the cross product with the desired axis. Cross product
  ! will be zero if the direction vector is the same as the axis

 !$AD II-LOOP
  do i=1,nCon
     resultDir = cross(axis,dirVec(i,:))

     x(i) = 0
     do j =1,3
        x(i) = x(i)+(resultDir(j))**2
     end do
     x(i) = sqrt(x(i))
  end do

end subroutine colinearity

subroutine planarity(axis, origin, p0, p1, p2, n, x, nCon)
  implicit none
  real, dimension(3)::origin, axis
  real, dimension(n,3)::p0,p1,p2
  real, dimension(nCon)::x
  integer:: n, nCon

  real, dimension(3*n,3)::allPoints,dist
  real, dimension(3*n)::scalarDist
  real :: tmp
  integer:: i, j

  ! copy data into all points array
  allPoints(1:n,:) = p0
  allPoints(n:2*n,:) = p1
  allPoints(2*n:3*n,:) = p2

  ! Compute the distance from the origin to each point
  x(1) = 0
  do i=1,n*3
     do j=1,3
        dist(i,j) = allPoints(i,j)-origin(j)
     end do

     !project it onto the axis
     !x(i) = dot_product(axis,dist(i,:))
     call dot(axis,dist(i,:),scalarDist(i))
     x(1) = x(1)+ scalarDist(i)**2!/(n*3)

  end do
  x(1) = sqrt(x(1))

end subroutine planarity

FUNCTION cross(a, b)
  real, DIMENSION(3) :: cross
  real, DIMENSION(3), INTENT(IN) :: a, b

  cross(1) = a(2) * b(3) - a(3) * b(2)
  cross(2) = a(3) * b(1) - a(1) * b(3)
  cross(3) = a(1) * b(2) - a(2) * b(1)
END FUNCTION cross

subroutine computeArea(p0, p1, p2, n, area)
    real, dimension (n,3):: p0, p1, p2, v1, v2
    real, dimension (n,3):: crosses
    real, dimension (n):: areas
    real, INTENT(OUT) :: area
    integer :: n
    v1 = 0.
    v2 = 0.
    crosses = 0.
    areas = 0.
    !$AD II-LOOP
    do i=1,n
        do j = 1,3
            v1(i,j) = p1(i,j) - p0(i,j)
            v2(i,j) = p2(i,j) - p0(i,j)
        end do
        crosses(i,:) = cross(v1(i,:),v2(i,:))
        do j = 1,3
            areas(i) = areas(i) + crosses(i,j)**2
        end do
        areas(i) = sqrt(areas(i))
    end do
    area = sum(areas)/2.
end subroutine

subroutine normalLengths(center, coords,n, X, nCon)

real, dimension(3)::center
real, dimension(n,3)::coords
real, dimension(nCon)::X
integer:: n,nCon

real ::refLength2,length2

refLength2 = 0
do i =1,3
 refLength2 = refLength2+(center(i)-coords(1,i))**2
end do

!$AD II-LOOP
do i =1,nCon
   length2 = 0
   do j =1,3
      Length2 = length2+(center(j)-coords(i+1,j))**2
   end do
   x(i) = sqrt(length2/reflength2)
end do
end subroutine normalLengths

SUBROUTINE PROJAREA(n, p0, p1, p2, axis, area)

  IMPLICIT NONE
! Input
  INTEGER, INTENT(IN) :: n
  REAL*8, INTENT(IN) :: p0(n, 3), p1(n, 3), p2(n, 3), axis(3)
! Output
  REAL*8, INTENT(OUT) :: area
! Working
  INTEGER :: i
  REAL*8 :: v1(3), v2(3), savec(3), pa
  DO i=1,n
    v1 = p1(i, :) - p0(i, :)
    v2 = p2(i, :) - p0(i, :)
    !CALL CROSS(v1, v2, savec)
    savec = CROSS(v1, v2)
    CALL DOT(savec, axis, pa)
    IF (pa .GT. 0) THEN
      area = area + pa
    ELSE
      pa = 0
    END IF
  END DO
  area = area/2.0

END SUBROUTINE PROJAREA

SUBROUTINE DOT(a, b, c)

  IMPLICIT NONE
  REAL*8, INTENT(IN) :: a(3), b(3)
  REAL*8, INTENT(OUT) :: c

  c = a(1)*b(1) + a(2)*b(2) + a(3)*b(3)

END SUBROUTINE DOT

max_stars_repo_name: camader/pygeo

max_issues_repo_name: camader/pygeo

max_forks_repo_name: camader/pygeo


text: While not raised in the city limits there are several local farms that raise Emus.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
        SUBROUTINE DSSKRC
      INCLUDE 'DSIOF.COM'
      INCLUDE 'XNSTRN.COM'
10      ID = IAND( IBASE( INDCLR ), MASKQ1 )
        IF ( ID .EQ. IDSRH ) GO TO 40
        IF ( ID .EQ. IDSSB ) GO TO 40
        IF ( ID .EQ. IDSEF ) GO TO 30
        IF ( ID .EQ. IDSEB ) GO TO 20
        CALL DSMSG ( 103 )
20      CALL DSRDNB
        GO TO 10
30      INDCLR = INDCLR + 1
        INDCBP = INDCLR
        GO TO 7000
40      LEN  = IAND( IBASE( INDCLR ), MASKH2 )
        ICLR = INDCLR + LEN + 1
        ID   = IAND( IBASE( ICLR ), MASKQ1 )
        IF ( ID .EQ. IDSRT ) GO TO 50
        CALL DSMSG ( 104 )
50      IFLG = IAND( IBASE( ICLR ), MASKQ2 )
        IF ( IFLG .EQ. IDSC ) GO TO 60
        CALL DSRDNB
        GO TO 10
60      INDCLR = ICLR +  1
        INDCBP = INDCLR
7000    RETURN
        END


!============================================================================
      SUBROUTINE grid_construction_left_low()
!============================================================================
! Region Left Low (LL)
!
! The specific internal energy along from Triple point to Critical
! point. The liquid region is bounded by the saturation line and the maximum
! pressure line, the minimum internal energy and the critical internal
! energy.
!
!  ** left isobar boundary (e,P)-->(T,v) 2D New_Raph
!     right saturation boundary (e)-->(vl,vv,T)-->(P) 3D New_Raph
!     middle (v,e)-->(T)--(P) 1D New_Raph
!
! The bilinear mapping is used for the interpolation from physical (v,e)
! doamin to the transformed domaine (X,Y).(Bicubic can be
! considered later by using also the derivatives)
!========================================================================
!
     USE def_constants
     USE def_variables
     USE non_linear_solvers
     USE properties
     USE grid_functions
!
      IMPLICIT NONE
!
      INTEGER  ::  i, j, N, exitflag, i1, i2, Niter
!      INTEGER, DIMENSION (NNN_LL,MMM_LL) ::  flag_pT_LL, N_pT_LL
!
      REAL(pr)  ::  u_min, res, delta, p_guess, extern, utest_max,&
&                   T_guess_sat,v_guess_sat, vap_guess_sat
!      REAL(pr)  :: p0, v_inters, T_inters, u_inters
!
      REAL(pr), DIMENSION (NNN_LL)  ::  vap_guess, vap_LL
!
      REAL(pr), DIMENSION (NNN_LL,MMM_LL)  ::  res_pT_LL, res_LL,&
      &                                        T_guess, v_guess,vvv_LL_log
      REAL(pr), DIMENSION (NNN_sat_LL) :: res_sat_LL,vap_sat_LL, T_sat_sat,&
&                                         T_sat_pmax,res_pmax_LL,T_liq_meta
      REAL(pr)  ::  energy,press, temp, out_2, sound,in_2,guess3,out3
!
! data used for bicubic interpolation
!
!REAL(8) :: v1, v2, v3, v4, u1, u2, u3, u4
!REAL(8) :: F1, F2, F3, F4, G1, G2, G3, G4
!REAL(8) :: a1, a2, a3, a4, b1, b2, b3, b4
!REAL(8) :: dF1v, dF2v, dF3v, dF4v, dF1u, dF2u, dF3u, dF4u, dF1vu, dF2vu,
!dF3vu, dF4vu
!REAL(8) :: dF1vv, dF2vv, dF3vv, dF4vv, dF1uu, dF2uu, dF3uu, dF4uu
!REAL(8) :: dF1x, dF2x, dF3x, dF4x, dF1y, dF2y, dF3y, dF4y, dF1xy, dF2xy,
!dF3xy, dF4xy
!
!REAL(8), DIMENSION (4, 4) :: A_bicub, C_bicub
!
!
! Initial guess for point (1,1) and (1,MMM_LL)
!
        T_guess(1,1) = T_tri
        v_guess(1,1) = 1_pr / rho_tri_L
!
        TTT_LL(1,MMM_LL) = T_tri
        T_guess(1,MMM_LL) = T_tri
!
        vvv_LL(1,MMM_LL) = 1_pr / rho_tri_L
        v_guess(1,MMM_LL) = 1_pr / rho_tri_L
        ppp_LL(1,MMM_LL) = P_tri
        vap_guess(1) = 1_pr / rho_tri_R
        vap_LL(1) = 1_pr / rho_tri_R
        ccc_LL(1,MMM_LL) = c_tri_L
!--------------------------------------------------------------------------
!
! Construct vector y_mesh_LL(physical domain), x_mesh_LL(tansformed
! domain)
!
!-------------------------------------------------------------------------
        utest_max = 1.0_pr*e_cr ! interal energy smaller than cirtical internal energy
        u_min = e_tri_L
!
        delta = (utest_max - u_min) / (NNN_LL - 1)
        y_mesh_LL = u_min + (/(i*delta,i=0,NNN_LL-1)/)
!        y_mesh_LL(NNN_LL) = utest_max
!
        delta = (x_mesh_max - x_mesh_min) / (MMM_LL -1)
        x_mesh_LL = x_mesh_min + (/(i*delta, i=0,MMM_LL-1)/)
!        x_mesh_LL(NNN_LL) = x_mesh_max
!
! array y_mesh_sat_LL
!
        delta             = (utest_max - u_min)/(NNN_sat_LL-1)
        y_mesh_sat_LL     = u_min + (/(i*delta, i=0,NNN_sat_LL-1)/)
!
!
!--------------------------------------------------------------------------
!
! Left boundary of LL domain. isobar Pmax,  and u_min-->u_max
!
!--------------------------------------------------------------------------
!
!P (iso) at each point of  left boundary
!
        ppp_LL(:,1) = p_max
!
        DO i=1,NNN_LL
!
              CALL New_Rap2D(2, TTT_LL(i,1), vvv_LL(i,1), &
              & res_LL(i,1), Niter, exitflag, ppp_LL(i,1), y_mesh_LL(i),&
              & T_guess(i,1), v_guess(i,1))
!
             CALL sound_speed(TTT_LL(i,1),vvv_LL(i,1), sound)
             ccc_LL(i,1) = sound
!
!             CALL inter_energy(TTT_LL(i,1),vvv_LL(i,1),energy)
!        print*, i, abs(y_mesh_LL(i) - energy ) / y_mesh_LL(i)
              IF (res_LL(i,1) > res_ref) THEN
                print*, "left_res IN LL", res_LL(i,1), "iter", Niter,"flag",exitflag
                STOP
              ENDIF
              T_guess(i+1,1) = TTT_LL(i,1)
              v_guess(i+1,1) = vvv_LL(i,1)
!
         END DO

!---------------------------------------------------------------------------
!
! Right boundary of LL domain. Saturation line u_min-->u_max
!
!--------------------------------------------------------------------------
        DO i=2, NNN_LL-1
              T_guess(i,MMM_LL) = TTT_LL(i-1,MMM_LL)
              v_guess(i,MMM_LL) = vvv_LL(i-1,MMM_LL)
              vap_guess(i) = vap_LL(i-1)
!
              CALL New_Rap3D(1,vvv_LL(i,MMM_LL),vap_LL(i),TTT_LL(i,MMM_LL),&
              & res_LL(i,MMM_LL), Niter, exitflag, y_mesh_LL(i),in_2,&
              & v_guess(i,MMM_LL), vap_guess(i),T_guess(i,MMM_LL))
!
              CALL pressure(TTT_LL(i,MMM_LL),vvv_LL(i,MMM_LL),press)
              ppp_LL(i,MMM_LL) = press
!
              CALL sound_speed(TTT_LL(i,MMM_LL),vvv_LL(i,MMM_LL), sound)
              ccc_LL(i,MMM_LL) = sound
!
!        print*, "vap_v", vap_LL(i)
              IF (res_LL(i,MMM_LL) > res_ref) THEN
                print*, "right_res IN LL", res_LL(i,MMM_LL), "iter", Niter
                print*, 'i',i,'v',vvv_LL(i,MMM_LL),'vap',vap_LL(i),'e',y_mesh_LL(i)
                STOP
              ENDIF
!
!        print*, "v_right", vvv_LL(i,MMM_LL), 't_right',TTT_LL(i,MMM_LL)
!        print*, 'e', y_mesh_LL()
        END DO
! Last point is fixed to critical point values
        vvv_LL(NNN_LL,MMM_LL) = 1_pr/rho_cr
        TTT_LL(NNN_LL,MMM_LL) = T_cr
        ppp_LL(NNN_LL,MMM_LL) = P_cr
        ccc_LL(NNN_LL,MMM_LL) = c_cr
!        print*, 'LL','n-1','e', y_mesh_LL(NNN_LL-1), vvv_LL(NNN_LL-1,MMM_LL),&
!&                               TTT_LL(NNN_LL-1,MMM_LL),ppp_LL(NNN_LL-1,MMM_LL),ccc_LL(NNN_LL-1,MMM_LL)
!        print*, 'LL','n',  'e', y_mesh_LL(NNN_LL), vvv_LL(NNN_LL,MMM_LL),&
!&                               TTT_LL(NNN_LL,MMM_LL),ppp_LL(NNN_LL,MMM_LL),ccc_LL(NNN_LL,MMM_LL)

!
!              T_guess(NNN_LL-2,MMM_LL) = TTT_LL(NNN_LL-3,MMM_LL)
!              v_guess(NNN_LL-2,MMM_LL) = 1_pr/rho_cr
!              vap_guess(NNN_LL-2) = vap_LL(NNN_LL-3)
!        DO i=NNN_LL-2, NNN_LL
!
!              CALL New_Rap3D(1,vvv_LL(i,MMM_LL),vap_LL(i),TTT_LL(i,MMM_LL),&
!              & res_LL(i,MMM_LL), Niter, exitflag, y_mesh_LL(i),&
!              & v_guess(i,MMM_LL), vap_guess(i),T_guess(i,MMM_LL))
!
!              CALL pressure(TTT_LL(i,MMM_LL),vvv_LL(i,MMM_LL),press)
!              ppp_LL(i,MMM_LL) = press
!
!             CALL sound_speed(TTT_LL(i,MMM_LL),vvv_LL(i,MMM_LL), sound)
!             ccc_LL(i,MMM_LL) = sound
!              T_guess(i+1,MMM_LL) = TTT_LL(i,MMM_LL)
!              v_guess(i+1,MMM_LL) = vvv_LL(i,MMM_LL)
!              vap_guess(i+1) = vap_LL(i)
!
!       print*, Niter, res_LL(i,MMM_LL)
!              IF (res_LL(i,MMM_LL) > res_ref) THEN
!                print*, "right_res", res_LL(i,MMM_LL), "iter", Niter
!                STOP
!              ENDIF
!
!        print*, "v_right", vvv_LL(i,MMM_LL), 't_right',TTT_LL(i,MMM_LL)
!        print*, 'e', y_mesh_LL(i)
!        END DO

!---------------------------------------------------------------------------
!
! Middle domain of the two boundaries. v_min-->v_max, u_min-->u_max
! vvv_LL in the physical space is built by the X in the transformed
! space with a LINEAR SCALING TRANSFORMATION
!
!--------------------------------------------------------------------------
        DO j=2,MMM_LL-1
           DO i=1, NNN_LL
!
                vvv_LL(i,j) = ((x_mesh_LL(j) - x_mesh_min) /&
                & (x_mesh_max - x_mesh_min)) *(vvv_LL(i,MMM_LL)-vvv_LL(i,1))+&
                & vvv_LL(i,1)
!                vvv_LL_log(i,j) = ((x_mesh_LL(j) - x_mesh_min) /&
!                & (x_mesh_max - x_mesh_min)) *(log10(vvv_LL(i,MMM_LL))-log10(vvv_LL(i,1)))+&
!                & log10(vvv_LL(i,1))
!               vvv_LL(i,j) = 10_pr**vvv_LL_log(i,j)
!
                CALL New_Rap1D(1, temp, out_2, res_LL(i,j), Niter,&
               & exitflag, y_mesh_LL(i), TTT_LL(i,j-1),&
               & vvv_LL(i,j), out3)
!
                CALL New_Rap1D(1, temp, out_2, res_LL(i,j), Niter,&
               & exitflag, y_mesh_LL(i), temp,&
               & vvv_LL(i,j), out3)

!
                TTT_LL(i,j) = temp
!
!
              CALL pressure(TTT_LL(i,j),vvv_LL(i,j),press)
              ppp_LL(i,j) = press
!
             CALL sound_speed(TTT_LL(i,j),vvv_LL(i,j), sound)
             ccc_LL(i,j) = sound
!
                IF (res_LL(i,j) > res_ref) THEN
                  print*, "middle_res_LL", res_LL(i,j), "iter", Niter,i,j
                  STOP
                ENDIF
!
!                CALL inter_energy(temp, vvv_LL(i,j),energy)
!                print*, "error e", (y_mesh_LL(i) - energy)/y_mesh_LL(i)
!                print*, "residu", res_LL(i,j)
           END DO
        END DO
!        print*, 'LL ','n-1 ','e ', y_mesh_LL(NNN_LL), vvv_LL(NNN_LL,MMM_LL-1),&
!&                                  TTT_LL(NNN_LL,MMM_LL-1),ppp_LL(NNN_LL,MMM_LL-1),ccc_LL(NNN_LL,MMM_LL-1)
!
!
!print*, 'grid Cp'
! grid for Cp
!
!        DO j = 1,MMM_LL
!           DO  i = 1, NNN_LL
!               CALL heat_cap_p(TTT_LL(i,j),vvv_LL(i,j),cpcp_LL(i,j))
!               IF ( (cpcp_LL(i,j) /= cpcp_LL(i,j)) .OR. (cpcp_LL(i,j) <= 0.0) .OR.&
!&                   (cpcp_LL(i,j) >  cp_cr  ) ) THEN
!!                  print*,'cp in LL kJ/(kgK)', cpcp_LL(i,j)/1000.0,i,j
!                  cpcp_LL(i,j) = cp_cr
!                ENDIF
!             ENDDO
!        ENDDO
!
!------------------------------------------------------------
!
!          MESH CONSTRUCTION for spline coefficients
!
!------------------------------------------------------------
!print*, 'start spline'
!
! SATURATION
!
        T_sat_sat(1) = T_tri
        v_Lsat_LL(1) = 1_pr/rho_tri_L
        vap_sat_LL(1)= 1_pr/rho_tri_R
! 21/08/2017 fix last point to critical point
        T_sat_sat(NNN_sat_LL) = T_cr
        v_Lsat_LL(NNN_sat_LL) = 1_pr/rho_cr
        vap_sat_LL(NNN_sat_LL)= 1_pr/rho_cr
!
        T_sat_sat(NNN_sat_LL-1) = 304.1281
        v_Lsat_LL(NNN_sat_LL-1) = 1_pr/rho_cr - 2d-5
        vap_sat_LL(NNN_sat_LL-1)= 1_pr/rho_cr + 2d-5
! 20/05/2018 fix (NNN_sat_LL-1) point value
        DO i=2, NNN_sat_LL-2
              T_guess_sat = T_sat_sat(i-1)+4d-6
              v_guess_sat = v_Lsat_LL(i-1)!+1d-6
              vap_guess_sat = vap_sat_LL(i-1)!-1d-6
!print*, i, T_guess_sat, v_guess_sat,vap_guess_sat
!
              CALL New_Rap3D(1,v_Lsat_LL(i),vap_sat_LL(i),T_sat_sat(i),&
              & res_sat_LL(i), Niter, exitflag, y_mesh_sat_LL(i),in_2,&
              & v_guess_sat, vap_guess_sat,T_guess_sat)

!              CALL New_Rap3D(1,v_Lsat_LL(i),vap_sat_LL(i),T_sat_sat(i),&
!              & res_sat_LL(i), Niter, exitflag, y_mesh_sat_LL(i),in_2,&
!              & v_Lsat_LL(i), vap_sat_LL(i), T_sat_sat(i))
!
!print*, i, T_sat_sat(i), v_Lsat_LL(i), vap_sat_LL(i)
!
              IF (res_sat_LL(i) > res_ref) THEN
                print*, "right_sat_res", res_sat_LL(i), "iter", Niter,&
              &         'T=',T_sat_sat(i),i
                print*, 'Tguess',T_guess_sat,'vlquess',v_guess_sat,vap_guess_sat
                print*,'IN LL CONSTRUCTION'
                STOP
              ENDIF
!
        END DO
!
!
! META: compute spinodal liquid between  -290;004 - -190.3 (kj/kg)
!
!e = -427.18 - -290.004  (kj/kg) p = 0
!      T_guess_sat = 263.0
!      v_guess_sat = 1.1e-3_pr
! DO i = 1, 172
!    CALL New_Rap2D(4, T_liq_meta(i), v_liq_meta(i), &
!              & res, N, exitflag, y_mesh_sat_LL(i),in_2,&
!              & T_guess_sat,v_guess_sat)
!!
!   IF (res > 1e-5) THEN
!       print*,'LL_sat_spin-Convergence failed calculating T_sat_sat at spinodal(',i,')'
!       print*,'resnorm     = ',res ,'N_iter = ',N
!       STOP
!   ENDIF
!      T_guess_sat = T_liq_meta(i)
!      v_guess_sat = v_liq_meta(i)+0.1e-4_pr
!      CALL pressure(T_liq_meta(i),v_liq_meta(i), press)
!   print*, 'meta T_liqmeta', T_liq_meta(i), 'vliqmeta', v_liq_meta(i), y_mesh_sat_LL(i),i
!   print*,' press', press
!    T_liq_meta(i) = T_sat_sat(i)
!    v_liq_meta(i) = v_Lsat_LL(i)
! ENDDO
!STOP
!
!
!
!
!e = -290.004 - -201.476 (kj/kg) dp/dv_T = 0
!      T_guess_sat = 275.0
!      v_guess_sat = v_tri_L +1e-4_pr
! DO i = 173, 284
!!
!        CALL New_Rap2D(3, T_liq_meta(i), v_liq_meta(i), &
!              & res, N, exitflag, y_mesh_sat_LL(i),in_2,&
!              & T_guess_sat,v_guess_sat)
!!
!   IF (res > 1e-5) THEN
!       print*,'LL_sat_spin-Convergence failed calculating T_sat_sat at spinodal(',i,')'
!       print*,'resnorm     = ',res ,'N_iter = ',N
!       STOP
!   ENDIF
!   IF (i<200) THEN
!      T_guess_sat = T_liq_meta(i)
!      v_guess_sat = v_liq_meta(i)+1e-4_pr
!!   ELSEIF (i>=284) THEN
!!      T_guess_sat = T_liq_meta(i)+0.1_pr
!!      v_guess_sat = v_liq_meta(i)
!   ELSE
!      T_guess_sat = T_liq_meta(i)
!      v_guess_sat = v_liq_meta(i)-1e-4_pr

!   ENDIF
!      CALL pressure(T_liq_meta(i),v_liq_meta(i), press)
!   print*, 'meta T_liqmeta', T_liq_meta(i), 'vliqmeta', v_liq_meta(i), y_mesh_sat_LL(i),i
!   print*,' press', press
!   print*, 'sat', v_Lsat_LL(i), T_sat_sat(i)
! ENDDO
!!e = -201.476 - -190.3 (kj/kg) saturation line
!   DO i = 285, NNN_sat_LL
!      T_liq_meta(i) = T_sat_sat(i)
!      v_liq_meta(i) = v_Lsat_LL(i)
! ENDDO
!!   v_liq_meta(NNN_sat_LL)  = 1_pr/rho_cr
!
!
!print*, 'start p_max spline'
! P_max
!
!
        T_guess(1,1) = T_tri
        v_guess(1,1) = 1_pr/rho_tri_L
!
        DO i=1,NNN_sat_LL
!
              CALL New_Rap2D(2, T_sat_pmax(i), v_Lpmax_LL(i), &
              & res_pmax_LL(i), Niter, exitflag, p_max, y_mesh_sat_LL(i),&
              & T_guess(i,1), v_guess(i,1))
!
              IF (res_pmax_LL(i) > res_ref) THEN
                print*, "left_pmax_res IN LL ", res_pmax_LL(i), "iter", Niter
                STOP
              ENDIF
              T_guess(i+1,1) = T_sat_pmax(i)
              v_guess(i+1,1) = v_Lpmax_LL(i)
!
         END DO

!
!
! construction of spline coefficients associated to saturation
! curve nodes(What I need), meta curve nodes and pmax curve nodes(What I  need)
!
DO i = 1, NNN_LL-1
   j = ord_spline * (i-1) + 1
   CALL polyfit(spline_Lsat_LL(:,i), y_mesh_sat_LL(j:j+ord_spline),v_Lsat_LL(j:j+ord_spline), ord_spline)
!   CALL polyfit(spline_Lsat_LL(:,i), y_mesh_sat_LL(j:j+3),LOG10(v_Lsat_LL(j:j+3)), ord_spline)
!   CALL polyfit(spline_meta_LL(:,i), y_mesh_sat_LL(j:j+ord_spline),v_liq_meta(j:j+ord_spline), ord_spline)
   CALL polyfit(spline_pmax_LL(:,i), y_mesh_sat_LL(j:j+ord_spline),v_Lpmax_LL(j:j+ord_spline), ord_spline)
ENDDO
!print*,'------------------------------------------------------'
print*,           'LL GRID FINISH'
!print*,'------------------------------------------------------'
!
!
!
      END SUBROUTINE grid_construction_left_low


text:
module arkode_interface
  use farkode_mod
  use farkode_arkstep_mod
  use fsunmatrix_dense_mod
  use fsunlinsol_dense_mod
  contains
  integer(c_int) function FARKDense(arkode_mem, N)
    use, intrinsic :: iso_c_binding
    use fnvector_serial_mod
    implicit none
    type(c_ptr),     value :: arkode_mem
    integer(c_long), value :: N
    type(c_ptr)            :: sunmat_A
    type(c_ptr)            :: sunlinsol_LS
    type(c_ptr)            :: sunvec_y
    integer(c_int)         :: ierr

    sunvec_y = FN_VNewEmpty_Serial(N)
    sunmat_A = FSUNDenseMatrix(N, N)
    sunlinsol_LS = FSUNDenseLinearSolver(sunvec_y, sunmat_A)
    ierr = FARKStepSetLinearSolver(arkode_mem, sunlinsol_LS, sunmat_A)
  end function FARKDense
end module arkode_interface

max_stars_repo_name: ghweber/amrex

max_issues_repo_name: ghweber/amrex

max_forks_repo_name: ghweber/amrex