text:
      subroutine qqb_gmgmjt_gs(p,msq)
************************************************************************
*    Author: J.M. Campbell                                             *
*    March, 2013.                                                      *
*    Matrix element SUBTRACTION squared averaged over initial colors   *
*    and spins                                                         *
c     q(-p1)+qbar(-p2) --> gam(p3) + gam(p4) + parton(p5) + parton(p6) *
************************************************************************

      implicit none
      include 'constants.f'
      include 'ptilde.f'
      include 'nflav.f'
      include 'qqgg.f'
      include 'phot_dip.f'
      include 'frag.f'
      include 'ewcharge.f'
      include 'msqbits.f'
      integer j,k,nd
c --- remember: nd will count the dipoles

      double precision p(mxpart,4),msq(maxd,-nf:nf,-nf:nf),facgg,
     & pdip(mxpart,4),pswap(mxpart,4)
      double precision
     & msq15_2(-nf:nf,-nf:nf),msq25_1(-nf:nf,-nf:nf),
     & msq16_2(-nf:nf,-nf:nf),msq26_1(-nf:nf,-nf:nf),
     & msq15_6(-nf:nf,-nf:nf),msq26_5(-nf:nf,-nf:nf),
     & msq16_5(-nf:nf,-nf:nf),msq25_6(-nf:nf,-nf:nf),
     & msq56_1v(-nf:nf,-nf:nf),msq56_2v(-nf:nf,-nf:nf),
     & msq26_5v(-nf:nf,-nf:nf),msq26_1v(-nf:nf,-nf:nf),
     & msq15_6v(-nf:nf,-nf:nf),msq16_2v(-nf:nf,-nf:nf),
     & msq16_5v(-nf:nf,-nf:nf),msq25_6v(-nf:nf,-nf:nf),
     & msq25_1v(-nf:nf,-nf:nf),
     & msq15_2v(-nf:nf,-nf:nf),
     & dummy(-nf:nf,-nf:nf),
     & sub15_2(4),sub25_1(4),sub16_2(4),sub26_1(4),
     & sub15_6(4),sub16_5(4),sub25_6(4),sub26_5(4),
     & sub56_1(4),sub56_2(4),sub56_1v,sub56_2v,
     & sub26_5v,sub25_1v,sub26_1v,sub16_5v,sub16_2v,sub15_2v,sub15_6v,
     & sub25_6v,
     & msq35_1(-nf:nf,-nf:nf),msq36_1(-nf:nf,-nf:nf),
     & msq45_1(-nf:nf,-nf:nf),msq46_1(-nf:nf,-nf:nf),
     & msq35_1_swap(-nf:nf,-nf:nf),msq36_1_swap(-nf:nf,-nf:nf),
     & msq45_1_swap(-nf:nf,-nf:nf),msq46_1_swap(-nf:nf,-nf:nf),
     & sub35_1,sub36_1,sub45_1,sub46_1,
     & msqbits35_1(12),msqbits36_1(12),msqbits45_1(12),msqbits46_1(12),
     & msqbits35_1_swap(12),msqbits36_1_swap(12),
     & msqbits45_1_swap(12),msqbits46_1_swap(12)
      external qqb_gmgmjt,qqb_gmgmjt_gvec,qqb_dirgam_g,qqb_dirgam_g_swap

      if (frag) then
        ndmax=10
      else
        ndmax=6
      endif

c--- calculate all the initial-initial dipoles
      call dips(1,p,1,5,2,sub15_2,sub15_2v,msq15_2,msq15_2v,
     . qqb_gmgmjt,qqb_gmgmjt_gvec)
      call dips(2,p,2,5,1,sub25_1,sub25_1v,msq25_1,msq25_1v,
     . qqb_gmgmjt,qqb_gmgmjt_gvec)
      call dips(3,p,1,6,2,sub16_2,sub16_2v,msq16_2,msq16_2v,
     . qqb_gmgmjt,qqb_gmgmjt_gvec)
      call dips(4,p,2,6,1,sub26_1,sub26_1v,msq26_1,msq26_1v,
     . qqb_gmgmjt,qqb_gmgmjt_gvec)

c--- now the basic initial final ones
      call dips(5,p,1,5,6,sub15_6,sub15_6v,msq15_6,msq15_6v,
     . qqb_gmgmjt,qqb_gmgmjt_gvec)
c--- called for final initial the routine only supplies new values for
c--- sub... and sub...v and msqv
      call dips(5,p,5,6,1,sub56_1,sub56_1v,dummy,msq56_1v,
     . qqb_gmgmjt,qqb_gmgmjt_gvec)
      call dips(5,p,1,6,5,sub16_5,sub16_5v,msq16_5,msq16_5v,
     . qqb_gmgmjt,qqb_gmgmjt_gvec)

      call dips(6,p,2,6,5,sub26_5,sub26_5v,msq26_5,msq26_5v,
     . qqb_gmgmjt,qqb_gmgmjt_gvec)
      call dips(6,p,5,6,2,sub56_2,sub56_2v,dummy,msq56_2v,
     . qqb_gmgmjt,qqb_gmgmjt_gvec)
      call dips(6,p,2,5,6,sub25_6,sub25_6v,msq25_6,msq25_6v,
     . qqb_gmgmjt,qqb_gmgmjt_gvec)
c--- extra photon dipoles in the case of fragmentation
c--- note: qqb_dirgam_g_swap is just a wrapper routine to
c---       qqb_dirgam_g, with p4 and p5 interchanged;
c---       it is used, for instance, to obtain contributions
c---       from q q~ -> q~ q instead of the usual ordering q q~ -> q q~
      if (frag) then
         msqbits(:)=zip
         call dipsfrag(7,p,3,5,1,sub35_1,msq35_1,qqb_dirgam_g)
         msqbits35_1(:)=msqbits(:)
         call fill_gmgmjt_swap(7,msq35_1_swap)
         msqbits35_1_swap(:)=msqbits(:)

         msqbits(:)=zip
         call dipsfrag(8,p,3,6,1,sub36_1,msq36_1,qqb_dirgam_g)
         msqbits36_1(:)=msqbits(:)
         call fill_gmgmjt_swap(8,msq36_1_swap)
         msqbits36_1_swap(:)=msqbits(:)

         msqbits(:)=zip
         call dipsfrag(9,p,4,5,1,sub45_1,msq45_1,qqb_dirgam_g)
         msqbits45_1(:)=msqbits(:)
         call fill_gmgmjt_swap(9,msq45_1_swap)
         msqbits45_1_swap(:)=msqbits(:)

         msqbits(:)=zip
         call dipsfrag(10,p,4,6,1,sub46_1,msq46_1,qqb_dirgam_g)
         msqbits46_1(:)=msqbits(:)
         call fill_gmgmjt_swap(10,msq46_1_swap)
         msqbits46_1_swap(:)=msqbits(:)

         phot_dip(7:10)=.true.
      endif

      do j=-nf,nf
      do k=-nf,nf
      do nd=1,ndmax
        msq(nd,j,k)=0d0
      enddo
      enddo
      enddo

c--- this is the factor to apply the right couplings
c--- in the gluon-gluon contribution
      facgg=
     &  (2d0*Q(2)**4+dfloat(nf-2)*Q(1)**4)
     & /(2d0*Q(2)**2+dfloat(nf-2)*Q(1)**2)

c--- 2-quark, 2-gluon subtraction pieces
      do j=-nf,nf
      do k=-nf,nf

      if ((j .ne. 0) .and. (k .ne. 0) .and. (j.ne.-k)) goto 19

c--- do only q-qb and qb-q cases
      if (  ((j .gt. 0).and.(k .lt. 0))
     . .or. ((j .lt. 0).and.(k .gt. 0))) then
C-----half=statistical factor
      msq(1,j,k)=-half*msq15_2(j,k)*sub15_2(qq)/xn
      msq(2,j,k)=-half*msq25_1(j,k)*sub25_1(qq)/xn
      msq(3,j,k)=-half*msq16_2(j,k)*sub16_2(qq)/xn
      msq(4,j,k)=-half*msq26_1(j,k)*sub26_1(qq)/xn
      msq(5,j,k)=half*xn*(
     .  msq15_6(j,k)*(sub15_6(qq)+0.5d0*sub56_1(gg))
     . +0.5d0*msq56_1v(j,k)*sub56_1v
     . +msq16_5(j,k)*(sub16_5(qq)+0.5d0*sub56_1(gg))
     . +0.5d0*msq56_1v(j,k)*sub56_1v)
      msq(6,j,k)=half*xn*(
     .  msq26_5(j,k)*(sub26_5(qq)+0.5d0*sub56_2(gg))
     . +0.5d0*msq56_2v(j,k)*sub56_2v
     . +msq25_6(j,k)*(sub25_6(qq)+0.5d0*sub56_2(gg))
     . +0.5d0*msq56_2v(j,k)*sub56_2v)
      elseif ((k .eq. 0).and.(j.ne.0)) then
c--- q-g and qb-g cases
      msq(2,j,k)=2d0*tr*msq25_1(j,-j)*sub25_1(qg)
      msq(3,j,k)=xn*msq16_2(j,k)*sub16_2(qq)
      msq(4,j,k)=xn*(msq26_1(j,k)*sub26_1(gg)+msq26_1v(j,k)*sub26_1v)
      msq(5,j,k)=-(msq16_5(j,k)*sub16_5(qq)+msq16_5(j,k)*sub56_1(qq))/xn
      msq(6,j,k)=xn*(msq26_5(j,k)*sub26_5(gg)+msq26_5v(j,k)*sub26_5v
     .              +msq26_5(j,k)*sub56_2(qq))
      if (frag) then
        msq(7,j,k)=Q(j)**2*msq35_1(j,k)*sub35_1/2d0
        msq(9,j,k)=Q(j)**2*msq45_1(j,k)*sub45_1/2d0
      endif

      elseif ((j .eq. 0).and.(k.ne.0)) then
c--- g-q and g-qb cases
      msq(1,j,k)=2d0*tr*msq15_2(-k,k)*sub15_2(qg)
      msq(3,j,k)=xn*(msq16_2(j,k)*sub16_2(gg)+msq16_2v(j,k)*sub16_2v)
      msq(4,j,k)=xn*msq26_1(j,k)*sub26_1(qq)
      msq(5,j,k)=xn*(msq16_5(j,k)*sub16_5(gg)+msq16_5v(j,k)*sub16_5v
     .              +msq15_6(j,k)*sub56_1(qq))
      msq(6,j,k)=-(msq26_5(j,k)*sub26_5(qq)+msq26_5(j,k)*sub56_2(qq))/xn
      if (frag) then
        msq(7,j,k)=Q(k)**2*msq35_1(j,k)*sub35_1/2d0
        msq(9,j,k)=Q(k)**2*msq45_1(j,k)*sub45_1/2d0
      endif

      elseif ((j .eq. 0).and.(k .eq. 0)) then
c--- g-g case (real process is g(p1)+g(p2) --> qb(p5)+q(p6)
c---Hence 15 split multiplies q(15)+g(p2)-->Z+q(p6)
c---Hence 25 split multiplies g(p1)+q(p25)-->Z+q(p6)
      msq(1,j,k)=(msq15_2(+1,k)+msq15_2(+2,k)+msq15_2(+3,k)
     .           +msq15_2(+4,k)+msq15_2(+5,k))*sub15_2(qg)*2d0*tr
      msq(2,j,k)=(msq25_1(k,+1)+msq25_1(k,+2)+msq25_1(k,+3)
     .           +msq25_1(k,+4)+msq25_1(k,+5))*sub25_1(qg)*2d0*tr
      msq(3,j,k)=(msq16_2(-5,k)+msq16_2(-4,k)+msq16_2(-3,k)
     .           +msq16_2(-2,k)+msq16_2(-1,k))*sub16_2(qg)*2d0*tr
      msq(4,j,k)=(msq26_1(k,-5)+msq26_1(k,-4)+msq26_1(k,-3)
     .           +msq26_1(k,-2)+msq26_1(k,-1))*sub26_1(qg)*2d0*tr
      if (frag) then
        msq(7,j,k)=msq35_1(j,k)*sub35_1/2d0*facgg
        msq(8,j,k)=msq36_1(j,k)*sub36_1/2d0*facgg
        msq(9,j,k)=msq45_1(j,k)*sub45_1/2d0*facgg
        msq(10,j,k)=msq46_1(j,k)*sub46_1/2d0*facgg
      endif

      endif

 19   continue
      enddo
      enddo

c--- 4-quark subtraction pieces
      do j=-nf,nf
      do k=-nf,nf

      if (((j .gt. 0).and.(k .gt. 0)) .or.
     .    ((j .lt. 0).and.(k .lt. 0))) then
c--q-q or qb-qb
      if (j.eq.k) then
      msq(1,j,k)=msq(1,j,k)+0.5d0*(xn-1d0/xn)
     .  *(msq15_2(0,k)*sub15_2(gq)+msq15_2v(0,k)*sub15_2v)
      msq(2,j,k)=msq(2,j,k)+0.5d0*(xn-1d0/xn)
     .  *(msq25_1(j,0)*sub25_1(gq)+msq25_1v(j,0)*sub25_1v)
      msq(3,j,k)=msq(3,j,k)+0.5d0*(xn-1d0/xn)
     .  *(msq16_2(0,k)*sub16_2(gq)+msq16_2v(0,k)*sub16_2v)
      msq(4,j,k)=msq(4,j,k)+0.5d0*(xn-1d0/xn)
     .  *(msq26_1(j,0)*sub26_1(gq)+msq26_1v(j,0)*sub26_1v)
      else
      msq(1,j,k)=msq(1,j,k)+(xn-1d0/xn)
     .  *(msq15_2(0,k)*sub15_2(gq)+msq15_2v(0,k)*sub15_2v)
      msq(4,j,k)=msq(4,j,k)+(xn-1d0/xn)
     .  *(msq26_1(j,0)*sub26_1(gq)+msq26_1v(j,0)*sub26_1v)
      endif
      if (frag) then
        msq(7,j,k)=Q(j)**2*msq35_1(j,k)*sub35_1/2d0
        msq(8,j,k)=Q(k)**2*msq36_1_swap(j,k)*sub36_1/2d0
        msq(9,j,k)=Q(j)**2*msq45_1(j,k)*sub45_1/2d0
        msq(10,j,k)=Q(k)**2*msq46_1_swap(j,k)*sub46_1/2d0
      endif
      elseif (((j .gt. 0).and.(k .lt. 0))
     &    .or.((j .lt. 0).and.(k .gt. 0)))  then
c q-qbar or qbar-q
      if (j.eq.-k) then
      msq(1,j,k)=msq(1,j,k)+(xn-1d0/xn)
     .  *(msq15_2(0,k)*sub15_2(gq)+msq15_2v(0,k)*sub15_2v)
      msq(4,j,k)=msq(4,j,k)+(xn-1d0/xn)
     .  *(msq26_1(j,0)*sub26_1(gq)+msq26_1v(j,0)*sub26_1v)
      msq(6,j,k)=msq(6,j,k)+2d0*tr*dfloat(nflav)
     .  *(msq26_5(j,k)*sub56_2(gq)-msq56_2v(j,k)*sub56_2v)
      if (frag) then
c--- note: subtraction terms use symmetry of qqb and qbq amplitudes
        if     ((abs(j) .eq. 2) .or. (abs(j) .eq. 4)) then
        msq(7,j,k)=sub35_1/2d0*(
     &   Q(2)**2*(msqbits35_1(uub_uub)+msqbits35_1(uub_ccb))
     &  +Q(1)**2*(3d0*msqbits35_1(uub_ddb)))
        msq(8,j,k)=sub36_1/2d0*(
     &   Q(2)**2*(msqbits36_1_swap(uub_uub)+msqbits36_1_swap(uub_ccb))
     &  +Q(1)**2*(3d0*msqbits36_1_swap(uub_ddb)))
        msq(9,j,k)=sub45_1/2d0*(
     &   Q(2)**2*(msqbits45_1(uub_uub)+msqbits45_1(uub_ccb))
     &  +Q(1)**2*(3d0*msqbits45_1(uub_ddb)))
        msq(10,j,k)=sub46_1/2d0*(
     &   Q(2)**2*(msqbits46_1_swap(uub_uub)+msqbits46_1_swap(uub_ccb))
     &  +Q(1)**2*(3d0*msqbits46_1_swap(uub_ddb)))
        elseif ((abs(j) .eq. 1) .or. (abs(j) .eq. 3)
     &     .or. (abs(j) .eq. 5)) then
        msq(7,j,k)=sub35_1/2d0*(
     &   Q(1)**2*(msqbits35_1(ddb_ddb)+2d0*msqbits35_1(ddb_ssb))
     &  +Q(2)**2*(2d0*msqbits35_1(ddb_uub)))
        msq(8,j,k)=sub36_1/2d0*(
     &   Q(1)**2*(msqbits36_1_swap(ddb_ddb)
     &           +2d0*msqbits36_1_swap(ddb_ssb))
     &  +Q(2)**2*(2d0*msqbits36_1_swap(ddb_uub)))
        msq(9,j,k)=sub45_1/2d0*(
     &   Q(1)**2*(msqbits45_1(ddb_ddb)+2d0*msqbits45_1(ddb_ssb))
     &  +Q(2)**2*(2d0*msqbits45_1(ddb_uub)))
        msq(10,j,k)=sub46_1/2d0*(
     &   Q(1)**2*(msqbits46_1_swap(ddb_ddb)
     &           +2d0*msqbits46_1_swap(ddb_ssb))
     &  +Q(2)**2*(2d0*msqbits46_1_swap(ddb_uub)))
        endif
      endif
      else
      msq(1,j,k)=msq(1,j,k)+(xn-1d0/xn)
     .  *(msq15_2(0,k)*sub15_2(gq)+msq15_2v(0,k)*sub15_2v)
      msq(4,j,k)=msq(4,j,k)+(xn-1d0/xn)
     .  *(msq26_1(j,0)*sub26_1(gq)+msq26_1v(j,0)*sub26_1v)
      if (frag) then
        if (j .gt. 0) then
        msq(7,j,k)=Q(j)**2*msq35_1(j,k)*sub35_1/2d0
        msq(8,j,k)=Q(k)**2*msq36_1_swap(j,k)*sub36_1/2d0
        msq(9,j,k)=Q(j)**2*msq45_1(j,k)*sub45_1/2d0
        msq(10,j,k)=Q(k)**2*msq46_1_swap(j,k)*sub46_1/2d0
        else
        msq(7,j,k)=Q(j)**2*msq35_1_swap(j,k)*sub35_1/2d0
        msq(8,j,k)=Q(k)**2*msq36_1(j,k)*sub36_1/2d0
        msq(9,j,k)=Q(j)**2*msq45_1_swap(j,k)*sub45_1/2d0
        msq(10,j,k)=Q(k)**2*msq46_1(j,k)*sub46_1/2d0
        endif
      endif
      endif
c--qbar-q
c      elseif ((j .lt. 0).and.(k .gt. 0)) then
c      if (j.eq.-k) then
c      msq(2,j,k)=msq(2,j,k)+(xn-1d0/xn)
c     .  *(msq25_1(j,0)*sub25_1(gq)+msq25_1v(j,0)*sub25_1v)
c      msq(3,j,k)=msq(3,j,k)+(xn-1d0/xn)
c     .  *(msq16_2(0,k)*sub16_2(gq)+msq16_2v(0,k)*sub16_2v)
c      msq(6,j,k)=msq(6,j,k)+2d0*tr*dfloat(nflav)
c     . *(msq26_5(j,k)*sub56_2(gq)-msq56_2v(j,k)*sub56_2v)
c      else
c      msq(2,j,k)=msq(2,j,k)+(xn-1d0/xn)
c     .  *(msq25_1(j,0)*sub25_1(gq)+msq25_1v(j,0)*sub25_1v)
c      msq(3,j,k)=msq(3,j,k)+(xn-1d0/xn)
c     .  *(msq16_2(0,k)*sub16_2(gq)+msq16_2v(0,k)*sub16_2v)
c      endif
      endif

      enddo
      enddo

      return
      end

      subroutine fill_gmgmjt_swap(nd,msq_swap)
c--- routine that calls qqb_dirgam_g_swap and fills matrix elements
c---  input: nd (dipole number)
c--- output:msq_swap (array of matrix elements after 4<->5 swap)
      implicit none
      include 'constants.f'
      include 'ptilde.f'
      include 'z_dip.f'
      include 'incldip.f'
      integer nd
      double precision pdip(mxpart,4),msq_swap(-nf:nf,-nf:nf)

      if (incldip(nd)) then
        call getptilde(nd,pdip)
        pdip(4,:)=pdip(4,:)/z_dip(nd)
        call qqb_dirgam_g_swap(pdip,msq_swap)
      else
        msq_swap(:,:)=0d0
      endif

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
Hi, Im still learning how to use the wiki.

20100926 08:55:04 nbsp Welcome to the wiki!  If you have questions, just ask.   Users/CovertProfessor

20100930 17:59:24 nbsp Hi Ian, and Welcome to the Wiki.  As I commented to Mega, Concierge looks like a good service.  When you have a moment, be sure to look over the Welcome to the Wiki/Business Owner Business Owner Welcome page.  Its got a lot of good information on the role of wiki pages about businesses, the tone, the informational (versus promotional) nature, and so on.  Good luck with the business! Users/TomGarberson

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      SUBROUTINE RDCMPS(LUN)

C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    RDCMPS
C   PRGMMR: WOOLLEN          ORG: NP20       DATE: 2000-09-19
C
C ABSTRACT: THIS SUBROUTINE UNCOMPRESSES AND UNPACKS THE NEXT SUBSET
C   FROM THE INTERNAL COMPRESSED MESSAGE BUFFER (ARRAY MBAY IN COMMON
C   BLOCK /BITBUF/) AND STORES THE UNPACKED SUBSET WITHIN THE INTERNAL
C   ARRAY VAL(*,LUN) IN COMMON BLOCK /USRINT/.
C
C PROGRAM HISTORY LOG:
C 2000-09-19  J. WOOLLEN -- ORIGINAL AUTHOR
C 2002-05-14  J. WOOLLEN -- IMPROVED GENERALITY, PREVIOUSLY RDCMPS
C                           WOULD NOT RECOGNIZE COMPRESSED DELAYED
C                           REPLICATION AS A LEGITIMATE DATA STRUCTURE
C 2003-11-04  J. ATOR    -- ADDED DOCUMENTATION
C 2003-11-04  S. BENDER  -- ADDED REMARKS/BUFRLIB ROUTINE
C                           INTERDEPENDENCIES
C 2003-11-04  D. KEYSER  -- MAXJL (MAXIMUM NUMBER OF JUMP/LINK ENTRIES)
C                           INCREASED FROM 15000 TO 16000 (WAS IN
C                           VERIFICATION VERSION); UNIFIED/PORTABLE FOR
C                           WRF; ADDED HISTORY DOCUMENTATION
C 2004-08-18  J. ATOR    -- INITIALIZE CVAL TO EMPTY BEFORE CALLING UPC;
C                           CORRECT LOGIC FOR WHEN A CHARACTER VALUE IS
C                           THE SAME FOR ALL SUBSETS IN A MESSAGE;
C                           MAXIMUM MESSAGE LENGTH INCREASED FROM
C                           20,000 TO 50,000 BYTES
C DART $Id$
C
C USAGE:    CALL RDCMPS (LUN)
C   INPUT ARGUMENT LIST:
C     LUN      - INTEGER: I/O STREAM INDEX INTO INTERNAL MEMORY ARRAYS
C
C REMARKS:
C    THIS ROUTINE CALLS:        UPB      UPC      USRTPL
C    THIS ROUTINE IS CALLED BY: READSB
C                               Normally not called by any application
C                               programs.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C   MACHINE:  PORTABLE TO ALL PLATFORMS
C
C$$$

      INCLUDE 'bufrlib.prm'

      COMMON /BITBUF/ MAXBYT,IBIT,IBAY(MXMSGLD4),MBYT(NFILES),
     .                MBAY(MXMSGLD4,NFILES)
      COMMON /MSGCWD/ NMSG(NFILES),NSUB(NFILES),MSUB(NFILES),
     .                INODE(NFILES),IDATE(NFILES)
      COMMON /TABLES/ MAXTAB,NTAB,TAG(MAXJL),TYP(MAXJL),KNT(MAXJL),
     .                JUMP(MAXJL),LINK(MAXJL),JMPB(MAXJL),
     .                IBT(MAXJL),IRF(MAXJL),ISC(MAXJL),
     .                ITP(MAXJL),VALI(MAXJL),KNTI(MAXJL),
     .                ISEQ(MAXJL,2),JSEQ(MAXJL)
      COMMON /USRINT/ NVAL(NFILES),INV(MAXJL,NFILES),VAL(MAXJL,NFILES)

      CHARACTER*10 TAG
      CHARACTER*8  CREF,CVAL
      CHARACTER*3  TYP
      EQUIVALENCE  (CVAL,RVAL)
      REAL*8       VAL,RVAL,UPS,TEN

      DATA TEN/10/

C-----------------------------------------------------------------------
C     Statement function to compute BUFR "missing value" for field
C     of length LBIT bits (all bits "on"):

      LPS(LBIT) = MAX(2**(LBIT)-1,1)

C     Statement function to decode the encoded BUFR value IVAL according
C     to the scale and reference values that are stored within index NODE
C     of the internal arrays ISC(*) and IRF(*):

      UPS(NODE) = (IVAL+IRF(NODE))*TEN**(-ISC(NODE))
C-----------------------------------------------------------------------

C  SETUP THE SUBSET TEMPLATE
C  -------------------------

      CALL USRTPL(LUN,1,1)

C  UNCOMPRESS A SUBSET INTO THE VAL ARRAY ACCORDING TO TABLE B
C  -----------------------------------------------------------

      NSBS = NSUB(LUN)

C     Note that we are going to unpack the (NSBS)th subset from within
C     the current BUFR message.

      IBIT = MBYT(LUN)

C     Loop through each element of the subset.

      N = 0

1     DO N=N+1,NVAL(LUN)
      NODE = INV(N,LUN)
      NBIT = IBT(NODE)
      ITYP = ITP(NODE)

C     In each of the following code blocks, the "local reference value"
C     for the element is determined first, followed by the 6-bit value
C     which indicates how many bits are used to store the increment
C     (i.e. offset) from this "local reference value".  Then, we jump
C     ahead to where this increment is stored for this particular subset,
C     unpack it, and add it to the "local reference value" to determine
C     the final uncompressed value for this element from this subset.

C     Note that, if an element has the same final uncompressed value
C     for each subset in the message, then the encoding rules for BUFR
C     compression dictate that the "local reference value" will be equal
C     to this value, the 6-bit increment length indicator will have
C     a value of zero, and the actual increments themselves will be
C     omitted from the message.

      IF(ITYP.EQ.1.OR.ITYP.EQ.2) THEN

C        This is a numeric element.

         CALL UPB(LREF,NBIT,MBAY(1,LUN),IBIT)
         CALL UPB(LINC,   6,MBAY(1,LUN),IBIT)
         JBIT = IBIT + LINC*(NSBS-1)
         CALL UPB(NINC,LINC,MBAY(1,LUN),JBIT)
         IF(NINC.EQ.LPS(LINC)) NINC = LPS(NBIT)
         IVAL = LREF+NINC
         IF(ITYP.EQ.1) THEN
            CALL USRTPL(LUN,N,IVAL)
            GOTO 1
         ENDIF
         IF(IVAL.LT.LPS(NBIT)) VAL(N,LUN) = UPS(NODE)
         IBIT = IBIT + LINC*MSUB(LUN)
      ELSEIF(ITYP.EQ.3) THEN

C        This is a character element.

         CREF = ' '
         CALL UPC(CREF,NBIT/8,MBAY(1,LUN),IBIT)
         CALL UPB(LINC,   6,MBAY(1,LUN),IBIT)
         IF(LINC.EQ.0) THEN
            CVAL = CREF
         ELSE
            JBIT = IBIT + LINC*(NSBS-1)*8
            CVAL = ' '
            CALL UPC(CVAL,LINC,MBAY(1,LUN),JBIT)
         ENDIF
         VAL(N,LUN) = RVAL
         IBIT = IBIT + 8*LINC*MSUB(LUN)
      ENDIF
      ENDDO

      RETURN
      END

max_stars_repo_name: hkershaw-brown/feature-preprocess

max_issues_repo_name: hkershaw-brown/feature-preprocess

max_forks_repo_name: hkershaw-brown/feature-preprocess


text:
	SUBROUTINE GR_CMPV  ( rmin, rmax, rint, maxlvl, nlvl, clvl,
     +			      iret )
C************************************************************************
C* GR_CMPV								*
C*									*
C* This subroutine defines contour levels, given the data range and	*
C* the contour interval.						*
C*									*
C* GR_CMPV  ( RMIN, RMAX, RINT, MAXLVL, NLVL, CLVL, IRET )		*
C*									*
C* Input parameters:							*
C*	RMIN		REAL		Minimum value			*
C*	RMAX		REAL		Maximum value			*
C*	RINT		REAL		Contour interval		*
C*	MAXLVL		INTEGER		Max number of contour levels	*
C*									*
C* Output parameters:							*
C*	NLVL		INTEGER		Number of contour levels	*
C*	CLVL (NLVL)	REAL		Contour levels			*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return 		*
C*					 -8 = invalid range		*
C**									*
C* Log:									*
C* M. desJardins/GSFC	 2/85						*
C* M. desJardins/GSFC	 9/88	Cleaned up for GEMPAK4			*
C* K. Brill/HPC		 8/02	Stop DO WHILE using KTMX		*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	REAL		clvl (*)
	LOGICAL		done
C*
	INCLUDE		'ERMISS.FNC'
C------------------------------------------------------------------------
	iret = 0
	nlvl = 0
C
C*	Check for errors in input.
C
	IF  ( ( rint .le. 0. ) .or. ( rmin .gt. rmax ) .or.
     +	      ( ERMISS ( rmin ) ) .or. ( ERMISS ( rmax ) ) .or.
     +	      ( ERMISS ( rint ) ) )  THEN
	    iret = -8
	    RETURN
	END IF
C
C*	Find multiplier to use for rint to compute minimum value.
C
	r   = ( rmin / rint )
	knt = INT ( r )
	IF  ( FLOAT (knt) .eq. knt )  knt = knt - 1
	IF  ( rmin .ge. 0. )  knt = knt + 1
	ktmx = INT ( rmax / rint ) + 1
C
C*	Add contour values until rmax is reached or until MAXLVL levels
C*	are generated.
C
	done = .false.
	DO WHILE  ( .not. done )
	    tval = knt * rint
	    IF  ( tval .ge. rmin )  THEN
		IF  ( tval .gt. rmax )  THEN
		    done = .true.
		  ELSE IF  ( nlvl .ge. maxlvl )  THEN
		    done = .true.
		  ELSE
		    nlvl = nlvl + 1
		    clvl ( nlvl ) = tval
		END IF
	    END IF
	    knt = knt + 1
	    IF ( knt .gt. ktmx ) done = .true.
	END DO
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
!
! Copyright (C) 2007 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
!----------------------------------------------------------------------------
SUBROUTINE g2_kin ( ik )
  !----------------------------------------------------------------------------
  !
  ! ... Calculation of kinetic energy - includes the case of the modified
  ! ... kinetic energy functional for variable-cell calculations
  !
  USE kinds,                ONLY : DP
  USE cell_base,            ONLY : tpiba2
  USE klist,                ONLY : xk
  USE gvect,                ONLY : g
  USE wvfct,                ONLY : g2kin, igk, npw, ecfixed, qcutz, q2sigma
  !
  IMPLICIT NONE
  !
  INTEGER, INTENT (IN) :: ik
  !
  ! ... local variables
  !
  INTEGER :: ig
  REAL(DP), EXTERNAL :: qe_erf
  !
  !
  g2kin(1:npw) = ( ( xk(1,ik) + g(1,igk(1:npw)) )**2 + &
                   ( xk(2,ik) + g(2,igk(1:npw)) )**2 + &
                   ( xk(3,ik) + g(3,igk(1:npw)) )**2 ) * tpiba2
  !
  IF ( qcutz > 0.D0 ) THEN
     !
     DO ig = 1, npw
        !
        g2kin(ig) = g2kin(ig) + qcutz * &
             ( 1.D0 + qe_erf( ( g2kin(ig) - ecfixed ) / q2sigma ) )
        !
     END DO
     !
  END IF
  !
  RETURN
  !
END SUBROUTINE g2_kin

max_stars_repo_name: utdsimmons/ohpc

max_issues_repo_name: utdsimmons/ohpc

max_forks_repo_name: utdsimmons/ohpc


text:
Address(Cascade Place) is a residential Culdesacs culdesac in East Davis.

 Intersecting Streets

  Tulip Lane


  subroutine preload_for_SYNCH
!
  implicit none
!
  integer :: nstep,nsurfmax,nlabel,ntheta,i
!
  double precision :: rmn,rmx,zmn,zmx,raxis,zaxis
  double precision, dimension(:),   allocatable :: rbeg,rsmall,qsaf,psisurf,phitor
  double precision, dimension(:,:), allocatable :: R_st,Z_st,bmod_st,sqgnorm_st
!
  open(1,file='preload_for_SYNCH.inp')
  read (1,*) nstep    !number of integration steps
  read (1,*) nlabel   !grid size over radial variabl
  read (1,*) ntheta   !grid size over poloidal angle
  read (1,*) nsurfmax !number of starting points between the
                      !magnetic axis and right box boundary
                      !when searching for the separatrix
  close(1)
!
  allocate(rbeg(nlabel),rsmall(nlabel),qsaf(nlabel),psisurf(nlabel),phitor(nlabel))
  allocate(R_st(nlabel,ntheta),Z_st(nlabel,ntheta),bmod_st(nlabel,ntheta),sqgnorm_st(nlabel,ntheta))
!
  call field_line_integration_for_SYNCH(nstep,nsurfmax,nlabel,ntheta,    &
                                        rmn,rmx,zmn,zmx,raxis,zaxis,     &
                                        rbeg,rsmall,qsaf,psisurf,phitor, &
                                        R_st,Z_st,bmod_st,sqgnorm_st)
!
  open(1,form='formatted',file='box_size_axis.dat')
  write (1,*) rmn,rmx, '<= rmn, rmx (cm)'
  write (1,*) zmn,zmx, '<= zmn, zmx (cm)'
  write (1,*) raxis,zaxis, '<= raxis, zaxis (cm)'
  close(1)
!
  open(1,form='formatted',file='flux_functions.dat')
  write (1,*) '# R_beg, r,  q, psi_pol, psi_tor'
! 2020-02-20: Lukas Bauer: rbeg is uninitialized and never used for computation. Still 'random' values are written in file.
  do i=1,nlabel
    write (1,*) rbeg(i),rsmall(i),qsaf(i),psisurf(i),phitor(i)
  enddo
  close(1)
!
  open(1,form='formatted',file='twodim_functions.dat')
  write (1,*) nlabel, ntheta, '<= nlabel, ntheta'
  write (1,*) 'R(label,theta)'
  do i=1,nlabel
    write (1,*) R_st(i,:)
  enddo
  write (1,*) 'Z(label,theta)'
  do i=1,nlabel
    write (1,*) Z_st(i,:)
  enddo
  write (1,*) 'B(label,theta)'
  do i=1,nlabel
    write (1,*) bmod_st(i,:)
  enddo
  write (1,*) 'sqrtg_norm(label,theta)'
  do i=1,nlabel
    write (1,*) sqgnorm_st(i,:)
  enddo
  close(1)
!
  deallocate(rbeg,rsmall,qsaf,psisurf,phitor)
  deallocate(R_st,Z_st,bmod_st,sqgnorm_st)
!
  end subroutine preload_for_SYNCH


      program slope_srtm30
! ===============================================
      implicit none

      integer                ::  ix, iy
      integer                ::  nx, ny
      parameter                 (nx=6000)  !! input
      parameter                 (ny=6000)

      real,allocatable       ::  elevtn(:,:)
      integer*1,allocatable  ::  water(:,:), osm(:,:)
      real,allocatable       ::  var(:,:)

      character*128          ::  rfile1, wfile1
      integer                ::  ios
      real                   ::  west1, south1
      character*7            ::  cname1

      character*128          ::  buf
! ===============================================
      allocate(elevtn(nx,ny),water(nx,ny),osm(nx,ny),var(nx,ny))

      call getarg(1,buf)
      read(buf,*) west1
      call getarg(2,buf)
      read(buf,*) south1
      call set_name(west1,south1,cname1)

      rfile1='./MERIT_DEM/'//trim(cname1)//'.bin'
      open(11,file=rfile1,form='unformatted',access='direct',recl=4*nx*ny,status='old',iostat=ios)
      if( ios/=0 ) stop
      read(11,rec=1) elevtn
      close(11)

      rfile1='./G3WBM/'//trim(cname1)//'_wat.bil'
      open(11,file=rfile1,form='unformatted',access='direct',recl=1*nx*ny)
      read(11,rec=1) water
      close(11)

      rfile1='../5deg/'//trim(cname1)//'.bil'
      open(11,file=rfile1,form='unformatted',access='direct',recl=1*nx*ny)
      read(11,rec=1) osm
      close(11)

! ===================================
      var(:,:)=0
      do iy=1, ny
        do ix=1, nx
          if( elevtn(ix,iy)==-9999 )then
            var(ix,iy)=99
          else
            if( water(ix,iy)>=50 )then
              if( osm(ix,iy)==0 ) var(ix,iy)=10
              if( osm(ix,iy)>=1 ) var(ix,iy)=11
            else
              if( osm(ix,iy)==0 ) var(ix,iy)=0
              if( osm(ix,iy)>=1 ) var(ix,iy)=1
            endif
          endif
        end do
      end do

      wfile1='./var_'//trim(cname1)//'.bin'
      open(21,file=wfile1,form='unformatted',access='direct',recl=4*nx*ny)
      write(21,rec=1) var
      close(21)

      end program slope_srtm30

      subroutine set_name(lon,lat,cname)
! ===============================================
      implicit none
!
      real            ::  lon, lat

      character*1     ::  ew, sn
      character*2     ::  clat
      character*3     ::  clon
      character*7     ::  cname
! ===============================================
      if( lon<0 )then
        ew='w'
        write(clon,'(i3.3)') int(-lon)
      else
        ew='e'
        write(clon,'(i3.3)')  int(lon)
      endif

      if( lat<0 )then
        sn='s'
        write(clat,'(i2.2)') int(-lat)
      else
        sn='n'
        write(clat,'(i2.2)')  int(lat)
      endif

      cname=sn//clat//ew//clon

      end subroutine set_name


text:
POINTER(a,b)
END PROGRAM

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf


text:
subroutine SHExpandGLQ(cilm, lmax, gridglq, w, plx, zero, norm, csphase, lmax_calc)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!	This program will expand a grid of data (regulary spaced in longitude,
!	and irregurarly in latitude according to the Guass-Legendre quadrature points)
!	into spherical harmonics by utilizing an FFT of each latitudinal band, and a
!	Guass-Legendre Quadrature in latitute. Note that the array PLX is
!	optional, and should not be precomputed when memory is an issue (i.e., LMAX>360).
!	It is implicitly assumed that the gridded data are bandlimited to degree LMAX.
!	If the optional parameter LMAX_CALC is specified, the spherical harmonic coefficients
!	will be calculated only up to and including this degree.
!
!	If PLX is not present, the Legendre functions are computed on the fly
!	using the scaling methodolgy presented in Holmes and Featherston (2002).
!	When NORM is 1,2 or 4, these are accurate to degree 2800. When NORM is 3, the
!	routine is only stable to about degree 15.
!
!	Calling Parameters:
!		IN
!			gridglq		Gridded data used in the expansion generated
!					from a call to MakeGridGLQ, with dimensions
!					(LMAX+1, 2*LMAX+1).
!			lmax		Spherical harmonic bandwidth of the input grid.
!					If LMAX_CALC is not specified, this corresponds
!					to the maximum spherical harmonic degree of the expansion.
!			w		Gauss-Legendre points used in the integrations
!					(determined from a call to PreCompute).
!		OUT
!			cilm 		Spherical harmonic coefficients of expansion with
!					dimensions (2, LMAX+1, LMAX+1), or if LMAX_CAL is
!					present (2, LMAX_CALC+1, LMAX_CALC+1).
!		OPTIONAL (IN)
!			plx		Input array of Associated Legendre Polnomials computed
!					at the Gauss-Legendre points (determined from a call to
!					PreCompute). If this is not included, then the optional
!					array ZERO MUST be inlcuded.
!			zero		Array of dimension LMAX+1 that contains the latitudinal
!					gridpoints used in the Gauss-Legendre quadrature integration
!					scheme, calculated from a call to PreCompute. This is only
!					needed if PLX is not given.
!			norm		Normalization to be used when calculating the Legendre functions
!						(1) "geodesy" (default)
!						(2) Schmidt
!						(3) unnormalized
!						(4) orthonormalized
!			csphase		1: Do not include the phase factor of (-1)^m
!					-1: Apply the phase factor of (-1)^m.
!			lmax_calc	The maximum spherical harmonic degree calculated in the
!					spherical harmonic expansion.
!
!	Dependencies:	FFTW3, NGLQSH, CSPHASE_DEFAULT
!
!	Written by Mark Wieczorek (October 2003)
!
!	September 3, 2005. 	Modifed so that the array plx is now optional.
!	September 26, 2005. 	Added optional argument NORM.
!	October 16, 2006. 	The Legendre functions are now computed within this program
!				during the summations over l and m. This leads to an increase in speed
!				of about a factor of 2.
!	July 23, 2007.		Added optional argument LMAX_CALC.
!	August 18, 2009.	Fixed bug that could have given incorrect results when LMAX_CAlC = 0
!	November 21, 2011	Fixed problem where saved variables used in Plm recursion were not recalculated
!				if NORM changed from one call to the next (with the same value of N).
!
!	Copyright (c) 2005-2011, Mark A. Wieczorek
!	All rights reserved.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	use FFTW3
	use SHTOOLS, only: NGLQSH, CSPHASE_DEFAULT

	implicit none
	real*8, intent(in) ::	w(:), gridglq(:,:)
	real*8, intent(in), optional ::	plx(:,:), zero(:)
	real*8, intent(out) ::	cilm(:,:,:)
	integer, intent(in) ::	lmax
	integer, intent(in), optional :: norm, csphase, lmax_calc
	integer :: 		nlong, nlat, i, l, m, k, kstart, l1, m1, phase, i_s, astat(4), &
				lnorm, lmax_comp
	real*8 :: 		pi, gridl(2*lmax+1), prod, scalef, rescalem, u, p, pmm, &
				pm1, pm2, z, fcoef1(2, lmax+1), fcoef2(2, lmax+1)
	complex*16 ::		cc(lmax+1)
	integer*8 :: 		plan
	real*8, save, allocatable ::	f1(:), f2(:), sqr(:), symsign(:)
	integer, save ::	lmax_old=0, norm_old = 0

	if (present(lmax_calc)) then
		lmax_comp = min(lmax, lmax_calc)
	else
		lmax_comp = lmax
	endif

	if (present(lmax_calc)) then
		if (lmax_calc > lmax) then
			print*, "Error --- SHExpandGLQ"
			print*, "LMAX_CALC must be less than or equal to the effective bandwidth of GRID."
			print*, "Effective bandwidth of GRID = N/2 - 1 = ", lmax
			print*, "LMAX_CALC = ", lmax_calc
			stop
		endif
	endif

	if (size(cilm(:,1,1)) < 2 .or. size(cilm(1,:,1)) < lmax_comp+1 .or. size(cilm(1,1,:)) < lmax_comp+1) then
		print*, "Error --- SHExpandGLQ"
		print*, "CILM must be dimensioned as (2, LMAX_COMP+1, LMAX_COMP+1) where"
		print*, "LMAX_COMP = MIN(LMAX+1, LMAX_CALC+1)"
		print*, "LMAX = ", lmax
		if (present(lmax_calc)) print*, "LMAX_CALC = ", lmax_calc
		print*, "Input dimension is ", size(cilm(:,1,1)), size(cilm(1,:,1)), size(cilm(1,1,:))
		stop
	elseif (size(gridglq(1,:)) < 2*lmax+1 .or. size(gridglq(:,1)) < lmax+1 ) then
		print*, "Error --- SHExpandGLQ"
		print*, "GRIDGLQ must be dimensioned as (LMAX+1, 2*LMAX+1) where LMAX is ", lmax
		print*, "Input array is dimensioned ", size(gridglq(:,1)), size(gridglq(1,:))
		stop
	elseif (size(w) < lmax+1) then
		print*, "Error --- SHExpandGLQ"
		print*, "W must be dimensioned as (LMAX+1) where LMAX is ", lmax
		print*, "Input array is dimensioned as ", size(w)
		stop
	endif

	if (present(plx)) then
		if (size(plx(:,1)) < lmax+1 .or. size(plx(1,:)) < (lmax+1)*(lmax+2)/2) then
			print*, "Error --- SHExpandGLQ"
			print*, "PLX must be dimensioned as (LMAX+1, (LMAX+1)*(LMAX+2)/2) where LMAX is ", lmax
			print*, "Input array is dimensioned as ", size(plx(:,1)), size(plx(1,:))
			stop
		endif
	elseif (present(zero)) then
		if (size(zero) < lmax+1) then
			print*, "Error --- SHExpandGLQ"
			print*, "ZERO must be dimensioned as (LMAX+1) where LMAX is ", lmax
			print*, "Input array is dimensioned ", size(zero)
			stop
		endif
	else
		print*, "Error --- SHExpandGLQ"
		print*, "Either PLX or ZERO must be specified."
		stop
	endif

	if (present(norm)) then
		if (norm > 4 .or. norm < 1) then
			print*, "Error - SHExpandGLQ"
			print*, "Parameter NORM must be 1 (geodesy), 2 (Schmidt), 3 (unnormalized), or 4 (orthonormalized)."
			print*, "Input value is ", norm
			stop
		endif

		lnorm = norm

	else
		lnorm = 1
	endif

	if (present(csphase)) then
     		if (csphase /= -1 .and. csphase /= 1) then
     			print*, "Error --- SHExpandGLQ"
     			print*, "CSPHASE must be 1 (exclude) or -1 (include)."
     			stop
     		else
     			phase = csphase
     		endif
     	else
     		phase = CSPHASE_DEFAULT
     	endif

	nlong = 2*(lmax+1) -1	! This is the number of points (and period) used
				! in the FFT in order that m=0 to lmax.

	nlat = NGLQSH(lmax)	! nlat = lmax+1

	pi = acos(-1.0d0)

	cilm = 0.0d0

	scalef = 1.0d-280

	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	!
	!	Calculate recursion constants used in computing the Legendre polynomials
	!
	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

	if ( (lmax_comp /= lmax_old .or. lnorm /= norm_old) .and. .not. present(plx) ) then

		if (allocated(sqr)) deallocate(sqr)
		if (allocated(f1)) deallocate(f1)
		if (allocated(f2)) deallocate(f2)
		if (allocated(symsign)) deallocate(symsign)

		allocate(sqr(2*lmax_comp+1), stat=astat(1))
		allocate(f1((lmax_comp+1)*(lmax_comp+2)/2), stat=astat(2))
		allocate(f2((lmax_comp+1)*(lmax_comp+2)/2), stat=astat(3))
		allocate(symsign((lmax_comp+1)*(lmax_comp+2)/2), stat=astat(4))

		if (astat(1) /= 0 .or. astat(2) /= 0 .or. astat(3) /= 0 .or. astat(4) /= 0) then
			print*, "Error --- SHExpandGLQ"
			print*, "Problem allocating arrays SQR, F1, F2, or SYMSIGN", astat(1), astat(2), astat(3), astat(4)
			stop
		endif

		!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
		!
		! 	Calculate signs used for symmetry of Legendre functions about equator
		!
		!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

		k = 0
		do l = 0, lmax_comp, 1
			do m = 0, l, 1
				k = k + 1
				if (mod(l-m,2) == 0) then
					symsign(k) = 1.0d0
				else
					symsign(k) = -1.0d0
				endif
			enddo
		enddo

		!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
		!
		!	Precompute square roots of integers that are used several times.
		!
		!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

		do l=1, 2*lmax_comp+1
			sqr(l) = sqrt(dble(l))
		enddo

		!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
		!
		! 	Precompute multiplicative factors used in recursion relationships
		! 		P(l,m) = x*f1(l,m)*P(l-1,m) - P(l-2,m)*f2(l,m)
		!		k = l*(l+1)/2 + m + 1
		!	Note that prefactors are not used for the case when m=l as a different
		!	recursion is used. Furthermore, for m=l-1, Plmbar(l-2,m) is assumed to be zero.
		!
		!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

		k = 1

		select case(lnorm)

			case(1,4)

				if (lmax_comp /= 0) then
					k = k + 1
					f1(k) = sqr(3)
					f2(k) = 0.0d0
					k = k + 1
				endif

				do l=2, lmax_comp, 1
					k = k + 1
					f1(k) = sqr(2*l-1) * sqr(2*l+1) / dble(l)
					f2(k) = dble(l-1) * sqr(2*l+1) / sqr(2*l-3) / dble(l)
					do m=1, l-2, 1
						k = k+1
						f1(k) = sqr(2*l+1) * sqr(2*l-1) / sqr(l+m) / sqr(l-m)
                				f2(k) = sqr(2*l+1) * sqr(l-m-1) * sqr(l+m-1) &
                  				 	/ sqr(2*l-3) / sqr(l+m) / sqr(l-m)
					enddo
					k = k+1
					f1(k) = sqr(2*l+1) * sqr(2*l-1) / sqr(l+m) / sqr(l-m)
                			f2(k) = 0.0d0
					k = k + 1
				enddo

			case(2)

				if (lmax_comp /= 0) then
					k = k + 1
					f1(k) = 1.0d0
					f2(k) = 0.0d0
					k = k + 1
				endif

				do l=2, lmax_comp, 1
					k = k + 1
					f1(k) = dble(2*l-1) /dble(l)
					f2(k) = dble(l-1) /dble(l)
					do m=1, l-2, 1
						k = k+1
						f1(k) = dble(2*l-1) / sqr(l+m) / sqr(l-m)
                  				f2(k) = sqr(l-m-1) * sqr(l+m-1) / sqr(l+m) / sqr(l-m)
					enddo
					k = k+1
					f1(k) = dble(2*l-1) / sqr(l+m) / sqr(l-m)
                  			f2(k) = 0.0d0
					k = k + 1
				enddo

			case(3)

				do l=1, lmax_comp, 1
					k = k + 1
					f1(k) = dble(2*l-1) /dble(l)
					f2(k) = dble(l-1) /dble(l)
					do m=1, l-1, 1
						k = k+1
						f1(k) = dble(2*l-1) / dble(l-m)
                  				f2(k) = dble(l+m-1) / dble(l-m)
					enddo
					k = k + 1
				enddo

		end select

		lmax_old = lmax_comp
		norm_old = lnorm

	endif

	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	!
	! 	Create generic plan for gridl
	!
	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

	call dfftw_plan_dft_r2c_1d_(plan, nlong, gridl(1:nlong), cc, FFTW_MEASURE)

	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	!
	!	Determine Cilms, one l at a time I by integrating over all
	!	latitudes using Gauss-Legendre Quadrature. When PLX is not
	!	present, the Legendre functions are computed on the fly
	!	during the summations over l and m. These are scaled using
	!	the methodology of Holmesand Featherstone (2002), with the
	!	exception of the m=0 terms that do not need to be scaled
	!
	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

	if (present(plx)) then

		do i=1, nlat

			gridl(1:nlong) = gridglq(i,1:nlong)
			call dfftw_execute_(plan)				! take fourier transform
			fcoef1(1,1:lmax+1) = 2.0d0 * w(i) * pi * dble(cc(1:lmax+1)) / dble(nlong)
			fcoef1(2,1:lmax+1) = -2.0d0 * w(i) * pi * dimag(cc(1:lmax+1)) / dble(nlong)

			k = 0
			do l = 0, lmax_comp, 1
				l1 = l + 1
				do m = 0, l, 1
					m1 = m+1
					k = k + 1
					cilm(1:2,l1,m1) = cilm(1:2,l1,m1) + plx(i,k) * fcoef1(1:2,m1)
				enddo
			enddo
		enddo

	else

		do i=1, (nlat+1)/2

			if (i == (nlat+1)/2 .and. mod(nlat,2) /= 0) then		! This latitude is the equator; z=0, u=1

				gridl(1:nlong) = gridglq(i,1:nlong)
				call dfftw_execute_(plan)				! take fourier transform
				fcoef1(1,1:lmax+1) = 2.0d0 * w(i) * pi * dble(cc(1:lmax+1)) / dble(nlong)
				fcoef1(2,1:lmax+1) = -2.0d0 * w(i) * pi * dimag(cc(1:lmax+1)) / dble(nlong)

				u = 1.0d0

				select case(lnorm)
					case(1,2,3);	pm2 = 1.0d0
					case(4);	pm2 = 1.0d0 / sqrt(4.0d0*pi)
				end select

				cilm(1,1,1) = cilm(1,1,1) + pm2 * fcoef1(1,1)

				if (lmax_comp == 0) cycle

				k = 2

				do l=2, lmax_comp, 2
					l1 = l+1
					k = k+l
					p = - f2(k) * pm2
					pm2 = p
					cilm(1,l1,1) = cilm(1,l1,1) + p * fcoef1(1,1)
					k = k + l + 1
				enddo

				select case(lnorm)
					case(1,2);	pmm = sqr(2) * scalef
					case(3);	pmm = scalef
					case(4);	pmm = sqr(2) * scalef / sqrt(4.0d0*pi)
				end select

				rescalem = 1.0d0/scalef
				kstart = 1

				do m = 1, lmax_comp-1, 1

					m1 = m + 1
					kstart = kstart+m+1

					select case(lnorm)
						case(1,4)
							pmm = phase * pmm * sqr(2*m+1) / sqr(2*m)
							pm2 = pmm
						case(2)
							pmm = phase * pmm * sqr(2*m+1) / sqr(2*m)
							pm2 = pmm / sqr(2*m+1)
						case(3)
							pmm = phase * pmm * dble(2*m-1)
							pm2 = pmm
					end select

					fcoef1(1:2,m1) = fcoef1(1:2,m1) * rescalem

					cilm(1:2,m1,m1) = cilm(1:2,m1,m1) + pm2 * fcoef1(1:2,m1)

					k = kstart+m+1

					do l = m+2, lmax_comp, 2
						l1 = l+1
						k = k+l
                  				p = - f2(k) * pm2
                  				pm2 = p
						cilm(1:2,l1,m1) = cilm(1:2,l1,m1) + p * fcoef1(1:2,m1)
						k = k + l + 1
					enddo

				enddo

            			select case(lnorm)
            				case(1,4);	pmm = phase * pmm * sqr(2*lmax_comp+1) / sqr(2*lmax_comp) * rescalem
            				case(2);	pmm = phase * pmm / sqr(2*lmax_comp) * rescalem
            				case(3);	pmm = phase * pmm * dble(2*lmax_comp-1) * rescalem
        			end select

        			cilm(1:2,lmax_comp+1,lmax_comp+1) = cilm(1:2,lmax_comp+1,lmax_comp+1) &
        					+ pmm * fcoef1(1:2,lmax_comp+1)

			else

				z = zero(i)
				u = sqrt( (1.0d0-z) * (1.0d0+z) )

				gridl(1:nlong) = gridglq(i,1:nlong)
				call dfftw_execute_(plan)
				fcoef1(1,1:lmax+1) = 2.0d0 * w(i) * pi * dble(cc(1:lmax+1)) / dble(nlong)
				fcoef1(2,1:lmax+1) = -2.0d0 * w(i) * pi * dimag(cc(1:lmax+1)) / dble(nlong)

				i_s = nlat + 1 - i	! point symmetric about the equator

				gridl(1:nlong) = gridglq(i_s,1:nlong)
				call dfftw_execute_(plan)
				fcoef2(1,1:lmax+1) = 2.0d0 * w(i_s) * pi * dble(cc(1:lmax+1)) / dble(nlong)
				fcoef2(2,1:lmax+1) = -2.0d0 * w(i_s) * pi * dimag(cc(1:lmax+1)) / dble(nlong)

				select case(lnorm)
					case(1,2,3);	pm2 = 1.0d0
					case(4);	pm2 = 1.0d0 / sqrt(4.0d0*pi)
				end select

				cilm(1,1,1) = cilm(1,1,1) + pm2 * (fcoef1(1,1) + fcoef2(1,1) )
				! symsign(1) = 1

				if (lmax_comp == 0) cycle

				k = 2
				pm1 =  f1(k) * z * pm2
				cilm(1,2,1) = cilm(1,2,1) + pm1 * ( fcoef1(1,1) - fcoef2(1,1) )
				! symsign(2) = -1

				do l=2, lmax_comp, 1
					l1 = l+1
					k = k+l
					p = f1(k) * z * pm1 - f2(k) * pm2
					pm2 = pm1
					pm1 = p
					cilm(1,l1,1) = cilm(1,l1,1) + p * ( fcoef1(1,1) + fcoef2(1,1) * symsign(k) )
				enddo

				select case(lnorm)
					case(1,2);	pmm = sqr(2) * scalef
					case(3);	pmm = scalef
					case(4);	pmm = sqr(2) * scalef / sqrt(4.0d0*pi)
				end select

				rescalem = 1.0d0/scalef
				kstart = 1

				do m = 1, lmax_comp-1, 1

					m1 = m+1
					rescalem = rescalem * u
					kstart = kstart+m+1

					select case(lnorm)
						case(1,4)
							pmm = phase * pmm * sqr(2*m+1) / sqr(2*m)
							pm2 = pmm
						case(2)
							pmm = phase * pmm * sqr(2*m+1) / sqr(2*m)
							pm2 = pmm / sqr(2*m+1)
						case(3)
							pmm = phase * pmm * dble(2*m-1)
							pm2 = pmm
					end select

					fcoef1(1:2,m1) = fcoef1(1:2,m1) * rescalem
					fcoef2(1:2,m1) = fcoef2(1:2,m1) * rescalem

					cilm(1:2,m1,m1) = cilm(1:2,m1,m1) + pm2 * &
							( fcoef1(1:2,m1) + fcoef2(1:2,m1) )
					! symsign(kstart) = 1

					k = kstart+m+1
	   				pm1 = z * f1(k) * pm2

	   				cilm(1:2,m1+1,m1) = cilm(1:2,m1+1,m1) + pm1 * &
               						( fcoef1(1:2,m1) - fcoef2(1:2,m1) )
               				! symsign(k) = -1

					do l = m+2, lmax_comp, 1
						l1 = l+1
						k = k+l
                  				p = z * f1(k) * pm1-f2(k) * pm2
                  				pm2 = pm1
                  				pm1 = p
						cilm(1:2,l1,m1) = cilm(1:2,l1,m1) + p * &
								( fcoef1(1:2,m1) + fcoef2(1:2,m1) * symsign(k) )
					enddo

				enddo

				rescalem = rescalem * u

            			select case(lnorm)
            				case(1,4);	pmm = phase * pmm * sqr(2*lmax_comp+1) / sqr(2*lmax_comp) * rescalem
            				case(2);	pmm = phase * pmm / sqr(2*lmax_comp) * rescalem
            				case(3);	pmm = phase * pmm * dble(2*lmax_comp-1) * rescalem
        			end select

        			cilm(1:2,lmax_comp+1,lmax_comp+1) = cilm(1:2,lmax_comp+1,lmax_comp+1) + pmm * &
        						( fcoef1(1:2,lmax_comp+1) + fcoef2(1:2,lmax_comp+1) )
        						! symsign = 1

			endif
		enddo

	endif

	call dfftw_destroy_plan_(plan)

	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	!
	! 	Divide by integral of Ylm*Ylm
	!
	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

	select case(lnorm)

		case(1)

			do l=0, lmax_comp, 1
				cilm(1:2,l+1, 1:l+1) = cilm(1:2,l+1, 1:l+1) / (4.0d0*pi)
			enddo

		case(2)

			do l=0, lmax_comp, 1
				cilm(1:2,l+1, 1:l+1) = dble(2*l+1) * cilm(1:2,l+1, 1:l+1) / (4.0d0*pi)
			enddo

		case(3)

			do l=0, lmax_comp, 1
				prod = 4.0d0*pi/dble(2*l+1)
				cilm(1,l+1,1) = cilm(1,l+1,1) / prod
				prod = prod / 2.0d0
				do m=1, l-1, 1
					prod = prod * dble(l+m) * dble(l-m+1)
					cilm(1:2,l+1,m+1) = cilm(1:2,l+1,m+1) / prod
				enddo
				!do m=l case
				if (l /= 0) cilm(1:2,l+1,l+1) = cilm(1:2,l+1, l+1)/(prod*dble(2*l))
			enddo

	end select

end subroutine SHExpandGLQ

max_stars_repo_name: JackWalpole/seismo-fortran-fork

max_issues_repo_name: JackWalpole/seismo-fortran-fork

max_forks_repo_name: JackWalpole/seismo-fortran-fork


text:
!-----------------------------------------------------------------------
!
! MODULE: plib_module
!> @brief
!> This module contains the variables that control parallel decomposition
!> and the subroutines for parallel environment setup. Only module that
!> requires MPI library interaction except for time_module (MPI_WTIME).
!
!-----------------------------------------------------------------------

MODULE plib_module

  USE global_module, ONLY: i_knd, r_knd, l_knd

  USE time_module, ONLY: wtime

  IMPLICIT NONE

  PUBLIC

  INTERFACE glmax
    MODULE PROCEDURE glmax_i, glmax_d, glmax_d_1d
  END INTERFACE glmax

  INTERFACE glmin
    MODULE PROCEDURE glmin_i, glmin_d
  END INTERFACE glmin

  INTERFACE bcast
    MODULE PROCEDURE bcast_i_scalar, bcast_i_1d, bcast_d_scalar,       &
      bcast_d_1d
  END INTERFACE bcast

  INTERFACE psend
    MODULE PROCEDURE psend_d_2d, psend_d_3d
  END INTERFACE psend

  INTERFACE precv
    MODULE PROCEDURE precv_d_2d, precv_d_3d
  END INTERFACE precv

  SAVE
!_______________________________________________________________________
!
! Module Input Variables
!
! npey     - Number of MPI processes in the y-direction
! npez     - Number of MPI processes in the z-direction
! ichunk   - Size of work chunks in the x-direction
! nthreads - Number of OpenMP threads for energy parallelism
! nnested  - number of nested threads
!_______________________________________________________________________

  INTEGER(i_knd) :: npey=1, npez=1, ichunk=4, nthreads=1, nnested=0
!_______________________________________________________________________
!
! Run-time variables
!
! Note: all ranks are zero based
!
! root     - root process for comm_snap, 0
!
! nproc    - Number of MPI processes
! iproc    - Rank of calling process in base communicator
!
! comm_snap  - base communicator, duplicated from MPI_COMM_WORLD
! comm_space - SDD communicator, ndimen-1 grid for 2-D (x-y) or
!              3-D (x-y-z) problems. Non-existent for 1-D (x) problems.
! sproc      - Rank of calling process in comm_space
!
! ycomm      - y-dimension process communicator
! zcomm      - z-dimension process communicator
! yproc      - PE column in SDD 2-D PE mesh (comm_space)
! zproc      - PE row in SDD 2-D PE mesh (comm_space)
! firsty     - logical determining if lowest yproc
! lasty      - logical determining if highest yproc
! firstz     - logical determining if lowest zproc
! lastz      - logical determining if highest zproc
! ylop       - rank of preceding yproc in ycomm
! yhip       - rank of succeeding yproc in ycomm
! zlop       - rank of preceding zproc in zcomm
! zhip       - rank of succeeding zproc in zcomm
!
! g_off      - group offset for message tags
!
! thread_level       - level of MPI thread support
! thread_single      - MPI_THREAD_SINGLE
! thread_funneled    - MPI_THREAD_FUNNELED
! thread_serialized  - MPI_THREAD_SERIALIZED
! thread_multiple    - MPI_THREAD_MULTIPLE
! lock               - OpenMP lock
!
! num_grth  - minimum number of nthreads and ng; used to ensure loop
!             over groups with communications is sized properly
! do_nested - true/false use nested threading, i.e., mini-KBA
!_______________________________________________________________________

  INTEGER(i_knd), PARAMETER :: root=0, g_off = 2**14

  INTEGER(i_knd) :: nproc, iproc, comm_snap, comm_space, sproc, ycomm, &
    zcomm, yproc, zproc, ylop, yhip, zlop, zhip, thread_level,         &
    thread_single, thread_funneled, thread_serialized, thread_multiple,&
    max_threads, num_grth

  LOGICAL(l_knd) :: firsty, lasty, firstz, lastz, do_nested

  INCLUDE 'mpif.h'

  INCLUDE 'omp_lib.h'

  INTEGER(omp_lock_kind) :: lock

  CONTAINS

  SUBROUTINE pinit ( t1 )

!-----------------------------------------------------------------------
!
! Initialize the MPI process environment. Replicate MPI_COMM_WORLD to
! local communicator. Get each process its rank and total size.
!
!-----------------------------------------------------------------------

    REAL(r_knd), INTENT(OUT) :: t1
!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: ierr
!_______________________________________________________________________
!
!   Set thread support levels
!_______________________________________________________________________

    thread_single     = MPI_THREAD_SINGLE
    thread_funneled   = MPI_THREAD_FUNNELED
    thread_serialized = MPI_THREAD_SERIALIZED
    thread_multiple   = MPI_THREAD_MULTIPLE
!_______________________________________________________________________
!
!   Initialize MPI and return available thread support level. Prefer
!   thread_multiple but require thread_serialized. Abort for insufficient
!   thread support. Don't want to use MPI_COMM_WORLD everywhere, so
!   duplicate to comm_snap. Start the timer.
!_______________________________________________________________________

    CALL MPI_INIT_THREAD ( thread_serialized, thread_level, ierr )

    CALL WTIME ( t1 )

    CALL MPI_COMM_DUP ( MPI_COMM_WORLD, comm_snap, ierr )
!_______________________________________________________________________
!
!   Get the communicator size and each process' rank within the main
!   communicator
!_______________________________________________________________________

    CALL MPI_COMM_SIZE ( comm_snap, nproc, ierr )
    CALL MPI_COMM_RANK ( comm_snap, iproc, ierr )
!_______________________________________________________________________
!
!   Put a barrier for every process to reach this point.
!_______________________________________________________________________

    CALL barrier ( comm_snap )
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE pinit

  SUBROUTINE barrier ( comm )

!-----------------------------------------------------------------------
!
! MPI barrier
!
!-----------------------------------------------------------------------

    INTEGER(i_knd), INTENT(IN) :: comm
!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: ierr
!_______________________________________________________________________

    CALL MPI_BARRIER ( comm, ierr )
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE barrier

  SUBROUTINE pcomm_set

!-----------------------------------------------------------------------
!
! Setup the SDD communicator.
!
!-----------------------------------------------------------------------
!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: ierr

    INTEGER(i_knd), DIMENSION(2) :: dims

    LOGICAL(l_knd) :: reorder

    LOGICAL(l_knd), DIMENSION(2) :: periodic, remain
!_______________________________________________________________________
!
!   Use MPI routines to create a Cartesian topology
!_______________________________________________________________________

    dims(1) = npez
    dims(2) = npey
    periodic = .FALSE.
    reorder = .TRUE.

    CALL MPI_CART_CREATE ( comm_snap, 2, dims, periodic, reorder,      &
      comm_space, ierr )
!_______________________________________________________________________
!
!   Set up the sub-communicators of the cartesian grid
!_______________________________________________________________________

    remain(1) = .FALSE.
    remain(2) = .TRUE.
    CALL MPI_CART_SUB ( comm_space, remain, ycomm, ierr )

    remain(1) = .TRUE.
    remain(2) = .FALSE.
    CALL MPI_CART_SUB ( comm_space, remain, zcomm, ierr )
!_______________________________________________________________________
!
!   Get comm_space, ycomm, and zcomm ranks
!_______________________________________________________________________

    CALL MPI_COMM_RANK ( comm_space, sproc, ierr )
    CALL MPI_COMM_RANK ( ycomm, yproc, ierr )
    CALL MPI_COMM_RANK ( zcomm, zproc, ierr )
!_______________________________________________________________________
!
!   Set some variables used during solution
!_______________________________________________________________________

    IF ( yproc == 0 ) THEN
      firsty = .TRUE.
      ylop = yproc
    ELSE
      firsty = .FALSE.
      ylop = yproc - 1
    END IF

    IF ( yproc == npey-1 ) THEN
      lasty = .TRUE.
      yhip = yproc
    ELSE
      lasty = .FALSE.
      yhip = yproc + 1
    END IF

    IF ( zproc == 0 ) THEN
      firstz = .TRUE.
      zlop = zproc
    ELSE
      firstz = .FALSE.
      zlop = zproc - 1
    END IF

    IF ( zproc == npez-1 ) THEN
      lastz = .TRUE.
      zhip = zproc
    ELSE
      lastz = .FALSE.
      zhip = zproc + 1
    END IF
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE pcomm_set

  SUBROUTINE pend

!-----------------------------------------------------------------------
!
! Call to end MPI processes
!
!-----------------------------------------------------------------------
!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: ierr
!_______________________________________________________________________

    CALL MPI_FINALIZE ( ierr )
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE pend

  SUBROUTINE glmax_i ( value, comm )

!-----------------------------------------------------------------------
!
! All reduce global max value (integer). Use specified communicator.
!
!-----------------------------------------------------------------------

    INTEGER(i_knd), INTENT(IN) :: comm

    INTEGER(i_knd), INTENT(INOUT) :: value
!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: ierr, x
!_______________________________________________________________________

    IF ( comm == MPI_COMM_NULL ) RETURN
    CALL MPI_ALLREDUCE ( value, x, 1, MPI_INTEGER, MPI_MAX, comm, ierr )
    value = x
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE glmax_i

  SUBROUTINE glmax_d ( value, comm )

!-----------------------------------------------------------------------
!
! All reduce global max value (double precision float). Use specified
! communicator.
!
!-----------------------------------------------------------------------

    INTEGER(i_knd), INTENT(IN) :: comm

    REAL(r_knd), INTENT(INOUT) :: value
!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: ierr

    REAL(r_knd) :: x
!_______________________________________________________________________

    IF ( comm == MPI_COMM_NULL ) RETURN
    CALL MPI_ALLREDUCE ( value, x, 1, MPI_DOUBLE_PRECISION, MPI_MAX,   &
      comm, ierr )
    value = x
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE glmax_d

  SUBROUTINE glmax_d_1d ( value, dlen, comm )

!-----------------------------------------------------------------------
!
! All reduce global max value (double precision float) for 1-d array.
! Use specified communicator.
!
!-----------------------------------------------------------------------

    INTEGER(i_knd), INTENT(IN) :: dlen, comm

    REAL(r_knd), DIMENSION(dlen), INTENT(INOUT) :: value
!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: ierr

    REAL(r_knd), DIMENSION(dlen) :: x
!_______________________________________________________________________

    IF ( comm == MPI_COMM_NULL ) RETURN
    CALL MPI_ALLREDUCE ( value, x, dlen, MPI_DOUBLE_PRECISION, MPI_MAX,&
      comm, ierr )
    value = x
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE glmax_d_1d

  SUBROUTINE glmin_i ( value, comm )

!-----------------------------------------------------------------------
!
! All reduce global min value (integer). Use specified communicator.
!
!-----------------------------------------------------------------------

    INTEGER(i_knd), INTENT(IN) :: comm

    INTEGER(i_knd), INTENT(INOUT) :: value
!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: ierr, x
!_______________________________________________________________________

    IF ( comm == MPI_COMM_NULL ) RETURN
    CALL MPI_ALLREDUCE ( value, x, 1, MPI_INTEGER, MPI_MIN, comm, ierr )
    value = x
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE glmin_i

  SUBROUTINE glmin_d ( value, comm )

!-----------------------------------------------------------------------
!
! All reduce global min value (double precision float). Use specified
! communicator.
!
!-----------------------------------------------------------------------

    INTEGER(i_knd), INTENT(IN) :: comm

    REAL(r_knd), INTENT(INOUT) :: value
!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: ierr

    REAL(r_knd) :: x
!_______________________________________________________________________

    IF ( comm == MPI_COMM_NULL ) RETURN
    CALL MPI_ALLREDUCE ( value, x, 1, MPI_DOUBLE_PRECISION, MPI_MIN,   &
      comm, ierr )
    value = x
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE glmin_d

  SUBROUTINE bcast_i_scalar ( value, comm, bproc )

!-----------------------------------------------------------------------
!
! Broadcast (integer scalar). Use specified communicator and casting
! proc.
!
!-----------------------------------------------------------------------

    INTEGER(i_knd), INTENT(IN) :: comm, bproc

    INTEGER(i_knd), INTENT(INOUT) :: value
!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: ierr
!_______________________________________________________________________

    IF ( nproc == 1 ) RETURN
    IF ( comm == MPI_COMM_NULL ) RETURN

    CALL MPI_BCAST ( value, 1, MPI_INTEGER, bproc, comm, ierr )
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE bcast_i_scalar

  SUBROUTINE bcast_i_1d ( value, ilen, comm, bproc )

!-----------------------------------------------------------------------
!
! Broadcast (integer 1d array). Use specified communicator and casting
! proc.
!
!-----------------------------------------------------------------------

    INTEGER(i_knd), INTENT(IN) :: ilen, comm, bproc

    INTEGER(i_knd), DIMENSION(ilen), INTENT(INOUT) :: value
!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: ierr
!_______________________________________________________________________

    IF ( nproc == 1 ) RETURN
    IF ( comm == MPI_COMM_NULL ) RETURN

    CALL MPI_BCAST ( value, ilen, MPI_INTEGER, bproc, comm,   &
      ierr )
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE bcast_i_1d

  SUBROUTINE bcast_d_scalar ( value, comm, bproc )

!-----------------------------------------------------------------------
!
! Broadcast (double precision float scalar). Use specified communicator
! and casting proc.
!
!-----------------------------------------------------------------------

    INTEGER(i_knd), INTENT(IN) :: comm, bproc

    REAL(r_knd), INTENT(INOUT) :: value
!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: ierr
!_______________________________________________________________________

    IF ( nproc == 1 ) RETURN
    IF ( comm == MPI_COMM_NULL ) RETURN

    CALL MPI_BCAST ( value, 1, MPI_DOUBLE_PRECISION, bproc, comm, ierr )
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE bcast_d_scalar

  SUBROUTINE bcast_d_1d ( value, dlen, comm, bproc )

!-----------------------------------------------------------------------
!
! Broadcast (double precision float 1d array). Use specified
! communicator and casting proc.
!
!-----------------------------------------------------------------------

    INTEGER(i_knd), INTENT(IN) :: dlen, comm, bproc

    REAL(r_knd), DIMENSION(dlen), INTENT(INOUT) :: value
!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: ierr
!_______________________________________________________________________

    IF ( nproc == 1 ) RETURN
    IF ( comm == MPI_COMM_NULL ) RETURN

    CALL MPI_BCAST ( value, dlen, MPI_DOUBLE_PRECISION, bproc, comm,   &
      ierr )
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE bcast_d_1d

  SUBROUTINE psend_d_2d ( proc, myproc, d1, d2, value, comm, mtag )

!-----------------------------------------------------------------------
!
! Send a rank-3 double precision array.
!
!-----------------------------------------------------------------------

    INTEGER(i_knd), INTENT(IN) :: proc, myproc, d1, d2, comm, mtag

    REAL(r_knd), DIMENSION(d1,d2), INTENT(IN) :: value
!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: dlen, ierr
!_______________________________________________________________________

    IF ( proc==myproc .OR. comm==MPI_COMM_NULL ) RETURN

    dlen = d1*d2

    CALL MPI_SEND ( value, dlen, MPI_DOUBLE_PRECISION, proc, mtag,     &
      comm, ierr )
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE psend_d_2d

  SUBROUTINE psend_d_3d ( proc, myproc, d1, d2, d3, value, comm, mtag )

!-----------------------------------------------------------------------
!
! Send a rank-3 double precision array.
!
!-----------------------------------------------------------------------

    INTEGER(i_knd), INTENT(IN) :: proc, myproc, d1, d2, d3, comm, mtag

    REAL(r_knd), DIMENSION(d1,d2,d3), INTENT(IN) :: value
!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: dlen, ierr
!_______________________________________________________________________

    IF ( proc==myproc .OR. comm==MPI_COMM_NULL ) RETURN

    dlen = d1*d2*d3

    CALL MPI_SEND ( value, dlen, MPI_DOUBLE_PRECISION, proc, mtag,     &
      comm, ierr )
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE psend_d_3d

  SUBROUTINE precv_d_2d ( proc, myproc, d1, d2, value, comm, mtag )

!-----------------------------------------------------------------------
!
! Send a rank-3 double precision array.
!
!-----------------------------------------------------------------------

    INTEGER(i_knd), INTENT(IN) :: proc, myproc, d1, d2, comm, mtag

    REAL(r_knd), DIMENSION(d1,d2), INTENT(INOUT) :: value
!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: dlen, ierr, istat(MPI_STATUS_SIZE)
!_______________________________________________________________________

    IF ( proc==myproc .OR. comm==MPI_COMM_NULL ) RETURN

    dlen = d1*d2

    CALL MPI_RECV ( value, dlen, MPI_DOUBLE_PRECISION, proc, mtag,     &
      comm, istat, ierr )
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE precv_d_2d

  SUBROUTINE precv_d_3d ( proc, myproc, d1, d2, d3, value, comm, mtag )

!-----------------------------------------------------------------------
!
! Send a rank-3 double precision array.
!
!-----------------------------------------------------------------------

    INTEGER(i_knd), INTENT(IN) :: proc, myproc, d1, d2, d3, comm, mtag

    REAL(r_knd), DIMENSION(d1,d2,d3), INTENT(INOUT) :: value
!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: dlen, ierr, istat(MPI_STATUS_SIZE)
!_______________________________________________________________________

    IF ( proc==myproc .OR. comm==MPI_COMM_NULL ) RETURN

    dlen = d1*d2*d3

    CALL MPI_RECV ( value, dlen, MPI_DOUBLE_PRECISION, proc, mtag,     &
      comm, istat, ierr )
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE precv_d_3d

  SUBROUTINE cartrank ( coord, rank, comm )

!-----------------------------------------------------------------------
!
! Return the rank of a proc defined by the coordinates of the Cartesian
! communicator.
!
!-----------------------------------------------------------------------

    INTEGER(i_knd), INTENT(OUT) :: rank

    INTEGER(i_knd), INTENT(IN) :: comm

    INTEGER(i_knd), DIMENSION(2), INTENT(IN) :: coord
!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: ierr
!_______________________________________________________________________

    IF ( comm == MPI_COMM_NULL ) RETURN

    CALL MPI_CART_RANK ( comm, coord, rank, ierr )
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE cartrank

  SUBROUTINE pinit_omp ( ierr, error )

!-----------------------------------------------------------------------
!
! Setup the number of OpenMP threads. Check if any proc is exceeding
! max threads. Reset and report if so.
!
!-----------------------------------------------------------------------

    CHARACTER(LEN=64), INTENT(OUT) :: error

    INTEGER(i_knd), INTENT(OUT) :: ierr
!_______________________________________________________________________
!
!   Local variables
!_______________________________________________________________________

    INTEGER(i_knd) :: max_threads
!_______________________________________________________________________

    ierr = 0
    error = ' '

    max_threads = OMP_GET_MAX_THREADS()

    IF ( nthreads > max_threads ) THEN
      ierr = 1
      nthreads = max_threads
    END IF

    CALL glmax ( ierr, comm_snap )
    IF ( ierr /= 0 ) THEN
      error = '*WARNING: PINIT_OMP: NTHREADS>MAX_THREADS; reset to' // &
              ' MAX_THREADS'
    END IF

    CALL OMP_SET_NUM_THREADS ( nthreads )
!_______________________________________________________________________
!
!   Setup for nested threading
!_______________________________________________________________________

    do_nested = nnested > 0
    CALL OMP_SET_NESTED ( do_nested )
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE pinit_omp

  SUBROUTINE plock_omp ( dowhat )

!-----------------------------------------------------------------------
!
! Operate on an OpenMP lock
!
!-----------------------------------------------------------------------

    CHARACTER(LEN=*), INTENT(IN) :: dowhat
!_______________________________________________________________________

    SELECT CASE ( dowhat )
      CASE ( 'init' )
        CALL OMP_INIT_LOCK ( lock )
      CASE ( 'set' )
        CALL OMP_SET_LOCK ( lock )
      CASE ( 'unset' )
        CALL OMP_UNSET_LOCK ( lock )
      CASE ( 'destroy' )
        CALL OMP_DESTROY_LOCK ( lock )
      CASE DEFAULT
        RETURN
    END SELECT
!_______________________________________________________________________
!_______________________________________________________________________

  END SUBROUTINE plock_omp

  FUNCTION thread_num ()

!-----------------------------------------------------------------------
!
! Return thread number of caller, [0, nthreads-1]. Maintains separation
! of main code and OpenMP by placing here.
!
!-----------------------------------------------------------------------

    INTEGER(i_knd) :: thread_num
!_______________________________________________________________________

    thread_num = OMP_GET_THREAD_NUM()
!_______________________________________________________________________
!_______________________________________________________________________

  END FUNCTION thread_num

END MODULE plib_module

max_stars_repo_name: arm-hpc/SNAP

max_issues_repo_name: arm-hpc/SNAP

max_forks_repo_name: arm-hpc/SNAP


text:
 T Shirts
New Bike Church Tshirts now available!

 Stickers


program test_matmul
    use utils
    use sub
    integer nn
    real(8), dimension(5000,5000) :: a, b, c
    real all_time1, all_time2
    real time_init
    real time21, time22
    real    ::  T1,T2
    real    ::  dt1, dt2, dt3
    all_time1 = sys_time()
    print *, 'Test matmul'
    call cpu_time(T1)
    CALL RANDOM_NUMBER(a)
    CALL RANDOM_NUMBER(b)
    call cpu_time(T2)
    dt1 = T2-T1
    print *, "init time = ", dt1
    print *, "         n","   cpu time", "         sys time", "      cpu/sys time"
    do nn = 500, 2500, 200
        time_init = sys_time()
        call cpu_time(T1)
        call matmul(a,b,c,nn)
        call cpu_time(T2)
        time_matmul = sys_time() - time_init
        dt2 = T2-T1

        print *, nn, dt2,  time_matmul, dt2/time_matmul
    end do
    all_time2 = sys_time()
    print *,"all time =", all_time2 - all_time1

end program test_matmul


text:
*DECK DSPIGMR
      SUBROUTINE DSPIGMR (NEQ, TN, Y, SAVF, B, WGHT, N, MAXL, MAXLP1,
     1  KMP, DELTA, HL0, JPRE, MNEWT, F, PSOL, NPSL, X, V, HES, Q,
     2  LGMR, WP, IWP, WK, DL, IFLAG)
      EXTERNAL F, PSOL
      INTEGER NEQ,N,MAXL,MAXLP1,KMP,JPRE,MNEWT,NPSL,LGMR,IWP,IFLAG
      DOUBLE PRECISION TN,Y,SAVF,B,WGHT,DELTA,HL0,X,V,HES,Q,WP,WK,DL
      DIMENSION NEQ(*), Y(*), SAVF(*), B(*), WGHT(*), X(*), V(N,*),
     1    HES(MAXLP1,*), Q(*), WP(*), IWP(*), WK(*), DL(*)
C-----------------------------------------------------------------------
C This routine solves the linear system A * x = b using a scaled
C preconditioned version of the Generalized Minimal Residual method.
C An initial guess of x = 0 is assumed.
C-----------------------------------------------------------------------
C
C      On entry
C
C          NEQ = problem size, passed to F and PSOL (NEQ(1) = N).
C
C           TN = current value of t.
C
C            Y = array containing current dependent variable vector.
C
C         SAVF = array containing current value of f(t,y).
C
C            B = the right hand side of the system A*x = b.
C                B is also used as work space when computing
C                the final approximation.
C                (B is the same as V(*,MAXL+1) in the call to DSPIGMR.)
C
C         WGHT = the vector of length N containing the nonzero
C                elements of the diagonal scaling matrix.
C
C            N = the order of the matrix A, and the lengths
C                of the vectors WGHT, B and X.
C
C         MAXL = the maximum allowable order of the matrix HES.
C
C       MAXLP1 = MAXL + 1, used for dynamic dimensioning of HES.
C
C          KMP = the number of previous vectors the new vector VNEW
C                must be made orthogonal to.  KMP .le. MAXL.
C
C        DELTA = tolerance on residuals b - A*x in weighted RMS-norm.
C
C          HL0 = current value of (step size h) * (coefficient l0).
C
C         JPRE = preconditioner type flag.
C
C        MNEWT = Newton iteration counter (.ge. 0).
C
C           WK = real work array used by routine DATV and PSOL.
C
C           DL = real work array used for calculation of the residual
C                norm RHO when the method is incomplete (KMP .lt. MAXL).
C                Not needed or referenced in complete case (KMP = MAXL).
C
C           WP = real work array used by preconditioner PSOL.
C
C          IWP = integer work array used by preconditioner PSOL.
C
C      On return
C
C         X    = the final computed approximation to the solution
C                of the system A*x = b.
C
C         LGMR = the number of iterations performed and
C                the current order of the upper Hessenberg
C                matrix HES.
C
C         NPSL = the number of calls to PSOL.
C
C         V    = the N by (LGMR+1) array containing the LGMR
C                orthogonal vectors V(*,1) to V(*,LGMR).
C
C         HES  = the upper triangular factor of the QR decomposition
C                of the (LGMR+1) by lgmr upper Hessenberg matrix whose
C                entries are the scaled inner-products of A*V(*,i)
C                and V(*,k).
C
C         Q    = real array of length 2*MAXL containing the components
C                of the Givens rotations used in the QR decomposition
C                of HES.  It is loaded in DHEQR and used in DHELS.
C
C        IFLAG = integer error flag:
C                0 means convergence in LGMR iterations, LGMR .le. MAXL.
C                1 means the convergence test did not pass in MAXL
C                  iterations, but the residual norm is .lt. 1,
C                  or .lt. norm(b) if MNEWT = 0, and so x is computed.
C                2 means the convergence test did not pass in MAXL
C                  iterations, residual .gt. 1, and X is undefined.
C                3 means there was a recoverable error in PSOL
C                  caused by the preconditioner being out of date.
C               -1 means there was a nonrecoverable error in PSOL.
C
C-----------------------------------------------------------------------
      INTEGER I, IER, INFO, IP1, I2, J, K, LL, LLP1
      DOUBLE PRECISION BNRM,BNRM0,C,DLNRM,PROD,RHO,S,SNORMW,DNRM2,TEM
C
      IFLAG = 0
      LGMR = 0
      NPSL = 0
C-----------------------------------------------------------------------
C The initial residual is the vector b.  Apply scaling to b, and test
C for an immediate return with X = 0 or X = b.
C-----------------------------------------------------------------------
      DO 10 I = 1,N
 10     V(I,1) = B(I)*WGHT(I)
      BNRM0 = DNRM2 (N, V, 1)
      BNRM = BNRM0
      IF (BNRM0 .GT. DELTA) GO TO 30
      IF (MNEWT .GT. 0) GO TO 20
      CALL DCOPY (N, B, 1, X, 1)
      RETURN
 20   DO 25 I = 1,N
 25     X(I) = 0.0D0
      RETURN
 30   CONTINUE
C Apply inverse of left preconditioner to vector b. --------------------
      IER = 0
      IF (JPRE .EQ. 0 .OR. JPRE .EQ. 2) GO TO 55
      CALL PSOL (NEQ, TN, Y, SAVF, WK, HL0, WP, IWP, B, 1, IER)
      NPSL = 1
      IF (IER .NE. 0) GO TO 300
C Calculate norm of scaled vector V(*,1) and normalize it. -------------
      DO 50 I = 1,N
 50     V(I,1) = B(I)*WGHT(I)
      BNRM = DNRM2 (N, V, 1)
      DELTA = DELTA*(BNRM/BNRM0)
 55   TEM = 1.0D0/BNRM
      CALL DSCAL (N, TEM, V(1,1), 1)
C Zero out the HES array. ----------------------------------------------
      DO 65 J = 1,MAXL
        DO 60 I = 1,MAXLP1
 60       HES(I,J) = 0.0D0
 65     CONTINUE
C-----------------------------------------------------------------------
C Main loop to compute the vectors V(*,2) to V(*,MAXL).
C The running product PROD is needed for the convergence test.
C-----------------------------------------------------------------------
      PROD = 1.0D0
      DO 90 LL = 1,MAXL
        LGMR = LL
C-----------------------------------------------------------------------
C Call routine DATV to compute VNEW = Abar*v(ll), where Abar is
C the matrix A with scaling and inverse preconditioner factors applied.
C Call routine DORTHOG to orthogonalize the new vector VNEW = V(*,LL+1).
C Call routine DHEQR to update the factors of HES.
C-----------------------------------------------------------------------
        CALL DATV (NEQ, Y, SAVF, V(1,LL), WGHT, X, F, PSOL, V(1,LL+1),
     1        WK, WP, IWP, HL0, JPRE, IER, NPSL)
        IF (IER .NE. 0) GO TO 300
        CALL DORTHOG (V(1,LL+1), V, HES, N, LL, MAXLP1, KMP, SNORMW)
        HES(LL+1,LL) = SNORMW
        CALL DHEQR (HES, MAXLP1, LL, Q, INFO, LL)
        IF (INFO .EQ. LL) GO TO 120
C-----------------------------------------------------------------------
C Update RHO, the estimate of the norm of the residual b - A*xl.
C If KMP .lt. MAXL, then the vectors V(*,1),...,V(*,LL+1) are not
C necessarily orthogonal for LL .gt. KMP.  The vector DL must then
C be computed, and its norm used in the calculation of RHO.
C-----------------------------------------------------------------------
        PROD = PROD*Q(2*LL)
        RHO = ABS(PROD*BNRM)
        IF ((LL.GT.KMP) .AND. (KMP.LT.MAXL)) THEN
          IF (LL .EQ. KMP+1) THEN
            CALL DCOPY (N, V(1,1), 1, DL, 1)
            DO 75 I = 1,KMP
              IP1 = I + 1
              I2 = I*2
              S = Q(I2)
              C = Q(I2-1)
              DO 70 K = 1,N
 70             DL(K) = S*DL(K) + C*V(K,IP1)
 75           CONTINUE
            ENDIF
          S = Q(2*LL)
          C = Q(2*LL-1)/SNORMW
          LLP1 = LL + 1
          DO 80 K = 1,N
 80         DL(K) = S*DL(K) + C*V(K,LLP1)
          DLNRM = DNRM2 (N, DL, 1)
          RHO = RHO*DLNRM
          ENDIF
C-----------------------------------------------------------------------
C Test for convergence.  If passed, compute approximation xl.
C if failed and LL .lt. MAXL, then continue iterating.
C-----------------------------------------------------------------------
        IF (RHO .LE. DELTA) GO TO 200
        IF (LL .EQ. MAXL) GO TO 100
C-----------------------------------------------------------------------
C Rescale so that the norm of V(1,LL+1) is one.
C-----------------------------------------------------------------------
        TEM = 1.0D0/SNORMW
        CALL DSCAL (N, TEM, V(1,LL+1), 1)
 90     CONTINUE
 100  CONTINUE
      IF (RHO .LE. 1.0D0) GO TO 150
      IF (RHO .LE. BNRM .AND. MNEWT .EQ. 0) GO TO 150
 120  CONTINUE
      IFLAG = 2
      RETURN
 150  IFLAG = 1
C-----------------------------------------------------------------------
C Compute the approximation xl to the solution.
C Since the vector X was used as work space, and the initial guess
C of the Newton correction is zero, X must be reset to zero.
C-----------------------------------------------------------------------
 200  CONTINUE
      LL = LGMR
      LLP1 = LL + 1
      DO 210 K = 1,LLP1
 210    B(K) = 0.0D0
      B(1) = BNRM
      CALL DHELS (HES, MAXLP1, LL, Q, B)
      DO 220 K = 1,N
 220    X(K) = 0.0D0
      DO 230 I = 1,LL
        CALL DAXPY (N, B(I), V(1,I), 1, X, 1)
 230    CONTINUE
      DO 240 I = 1,N
 240    X(I) = X(I)/WGHT(I)
      IF (JPRE .LE. 1) RETURN
      CALL PSOL (NEQ, TN, Y, SAVF, WK, HL0, WP, IWP, X, 2, IER)
      NPSL = NPSL + 1
      IF (IER .NE. 0) GO TO 300
      RETURN
C-----------------------------------------------------------------------
C This block handles error returns forced by routine PSOL.
C-----------------------------------------------------------------------
 300  CONTINUE
      IF (IER .LT. 0) IFLAG = -1
      IF (IER .GT. 0) IFLAG = 3
C
      RETURN
C----------------------- End of Subroutine DSPIGMR ---------------------
      END

max_stars_repo_name: sbwhitecap/clocc-hg

max_issues_repo_name: sbwhitecap/clocc-hg

max_forks_repo_name: sbwhitecap/clocc-hg


text:
      SUBROUTINE E02054 (USEPRM)

C  E02054 tests the handling of error number 54

      COMMON /GLOBNU/ CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT, DUMRL
      INTEGER         CTLHND, ERRSIG, ERRFIL, IERRCT, UNERR,
     1        TESTCT, IFLERR, PASSSW, ERRSW, MAXLIN,
     2        CONID, MEMUN, WKID, WTYPE, GLBLUN, INDLUN,
     3        DUMINT(20), ERRIND
      REAL    DUMRL(20)

      COMMON /ERRINF/ ERRCOM,FUNCOM,FILCOM, ERNMSW, EXPSIZ,EXPERR,
     1                USRERR, ERRSAV,      FUNSAV,      FILSAV,
     2                EFCNT, EFID
      INTEGER         ERRCOM,FUNCOM,FILCOM, ERNMSW, EXPSIZ,EXPERR(10),
     1                USRERR, ERRSAV(200), FUNSAV(200), FILSAV(200),
     2                EFCNT, EFID(100)
      COMMON /ERRCHR/ CURCON,     ERRSRS,    ERRMRK,    ERFLNM,
     1                CONTAB
      CHARACTER       CURCON*200, ERRSRS*40, ERRMRK*20, ERFLNM*80,
     1                CONTAB(40)*150

C control flag
C                conditionally  always
      INTEGER    PCONDI,        PALWAY
      PARAMETER (PCONDI = 0,    PALWAY = 1)
C regeneration flag
      INTEGER    PPOSTP,     PPERFO
      PARAMETER (PPOSTP = 0, PPERFO = 1)
C text precision
      INTEGER    PSTRP,     PCHARP,     PSTRKP
      PARAMETER (PSTRP = 0, PCHARP = 1, PSTRKP = 2)
C clipping indicator
C                noclip      clip
      INTEGER    PNCLIP,     PCLIP
      PARAMETER (PNCLIP = 0, PCLIP = 1)
C modification mode
      INTEGER    PNIVE,     PUWOR,     PUQUM
      PARAMETER (PNIVE = 0, PUWOR = 1, PUQUM = 2)
C relative input priority
      INTEGER     PHIGHR,     PLOWER
      PARAMETER  (PHIGHR=0,   PLOWER=1)
C deferring model
      INTEGER     PASAP,     PBNIG,    PBNIL,    PAST1,    PWAIT
      PARAMETER  (PASAP=0,   PBNIG=1,  PBNIL=2,  PAST1=3,  PWAIT=4)

      INTEGER   USEPRM, COLIA(1,1), ISIZ, XYCLPI, BCLIPI
      INTEGER   FCLIPI, NUMOWK, OPWKID
      REAL      VWORMT(4,4), VWMPMT(4,4), VWCPLM(6), CSPEC(3)
      LOGICAL   STREQ

      CURCON = 'the specified workstation is not open'

      CALL SETVS ('54,3', EXPERR, EXPSIZ)
      ERRSRS = '4'
      CALL ESETUP (USEPRM)

C   generate errors
      CALL RCLWK (WKID)
      CALL TSTIGN (STREQ('OC**'))

      CALL RRST (WKID, PALWAY)
      CALL TSTIGN (STREQ('OC**'))

      CALL RUWK (WKID, PPERFO)
      CALL TSTIGN (STREQ('OC**'))

      CALL RSTXR (WKID, 1, 1, PSTRP, 1., 0.1, 1)
      CALL TSTIGN (STREQ('OC**'))

      CALL SETVS ('1', COLIA, ISIZ)
      CALL RSPAR (WKID, 1,1,1,1,1,1,1, COLIA)

      CALL SETRVS ('1.,0.,0.', CSPEC, ISIZ)
      CALL RSCR (WKID, 1, 3, CSPEC)

      CALL RSHRM (WKID, 0)

      CALL IDMAT (4, VWORMT)
      CALL IDMAT (4, VWMPMT)
      CALL SETRVS ('0.,1.,0.,1.,0.,1.', VWCPLM, ISIZ)
      XYCLPI = PCLIP
      BCLIPI = PCLIP
      FCLIPI = PCLIP
      CALL RSVWR3 (WKID, 1, VWORMT, VWMPMT, VWCPLM, XYCLPI, BCLIPI,
     1             FCLIPI)

      CALL RSVTIP (WKID, 1, 1, PHIGHR)

      CALL POPST (101)
      CALL RUPOST (WKID, 101)
      CALL TSTIGN (STREQ('OCO*'))
      CALL PCLST

      CALL POPWK (WKID, CONID, WTYPE)

      CALL RRST (WKID+1, PALWAY)
      CALL PQOPWK (1, ERRIND, NUMOWK, OPWKID)
      CALL CHKINQ ('pqopwk', ERRIND)
      CALL TSTIGN (STREQ('OO**') .AND. NUMOWK.EQ.1 .AND.
     1             OPWKID.EQ.WKID)

      CALL RSDUS (WKID+1, PASAP, PNIVE)

      CALL RSPAR (WKID+1, 1,1,1,1,1,1,1, COLIA)

      CALL RSWKW (WKID+1, 0.,1.,0.,1.)
      CALL PQOPWK (1, ERRIND, NUMOWK, OPWKID)
      CALL CHKINQ ('pqopwk', ERRIND)
      CALL TSTIGN (STREQ('OO**') .AND. NUMOWK.EQ.1 .AND.
     1             OPWKID.EQ.WKID)

      CALL POPST (101)
      CALL RPOST (WKID+1, 101, 1.0)
      CALL PQOPWK (1, ERRIND, NUMOWK, OPWKID)
      CALL CHKINQ ('pqopwk', ERRIND)
      CALL TSTIGN (STREQ('OOO*') .AND. NUMOWK.EQ.1 .AND.
     1             OPWKID.EQ.WKID)

      CALL PCLST
      CALL PCLWK (WKID)

C  examine saved errors

      CALL ENDERR

      END


text:
module init_data_module

  use multifab_module

  implicit none

  double precision, parameter :: Pi = 3.141592653589793238462643383279502884197d0

  private

  public :: init_data

contains

  subroutine init_data(data,dx,plo,phi)

    type(multifab),   intent(inout) :: data
    double precision, intent(in   ) :: dx(data%dim)
    double precision, intent(in   ) :: plo(data%dim), phi(data%dim)

    integer                   :: lo(data%dim), hi(data%dim), ng, i
    double precision, pointer :: dp(:,:,:,:)

    ng = data%ng

    do i=1,nfabs(data)
       dp => dataptr(data,i)
       lo = lwb(get_box(data,i))
       hi = upb(get_box(data,i))

       select case(data%dim)
       case (2)
          call bl_error('We only support 3-D')
       case (3)
          call init_data_3d(lo,hi,ng,dx,dp,plo,phi)
       end select
    end do

  end subroutine init_data

  subroutine init_data_3d(lo,hi,ng,dx,cons,phlo,phhi)

    use variables_module, only : irho, imx,imy,imz,iene,iry1,ncons, iH2, iO2, iN2
    use chemistry_module, only : nspecies, patm
    use probin_module,    only : rfire

    integer,          intent(in   ) :: lo(3),hi(3),ng
    double precision, intent(in   ) :: dx(3),phlo(3),phhi(3)
    double precision, intent(inout) :: cons(-ng+lo(1):hi(1)+ng,-ng+lo(2):hi(2)+ng,-ng+lo(3):hi(3)+ng,ncons)

    integer          :: i,j,k,n
    double precision :: x, y, z, r

    double precision Xt(nspecies), Yt(nspecies)
    double precision rhot,u1t,u2t,u3t,Tt,et,Pt
    integer :: iwrk
    double precision :: rwrk, expfac, kx, ky, kz

    kx = 2.d0*Pi/(phhi(1) - phlo(1))
    ky = 2.d0*Pi/(phhi(2) - phlo(2))
    kz = 2.d0*Pi/(phhi(3) - phlo(3))

    !$omp parallel do &
    !$omp private(i,j,k,n,x,y,z,r) &
    !$omp private(Xt,Yt,rhot,u1t,u2t,u3t,Tt,et,Pt,iwrk,rwrk,expfac)
    do k=lo(3),hi(3)
       z = phlo(3) + dx(3)*(k + 0.5d0)
       do j=lo(2),hi(2)
          y = phlo(2) + dx(2)*(j + 0.5d0)
          do i=lo(1),hi(1)
             x = phlo(1) + dx(1)*(i + 0.5d0)

             r = sqrt(x**2+y**2+z**2)

             Pt = Patm
             Tt = 300.0d0

             Xt = 0.0d0
             Xt(iH2) = 0.10d0
             Xt(iO2) = 0.25d0

             expfac = exp(-(r / rfire)**2)
             Pt      = Pt      + 0.1d0*patm * expfac
             Tt      = Tt      + 1100.0d0   * expfac
             Xt(iH2) = Xt(iH2) + 0.025d0    * expfac
             Xt(iO2) = Xt(iO2) - 0.050d0    * expfac

             Xt(iN2) = 1.0d0 - Xt(iH2) - Xt(iO2)

             u1t =  sin(kx*x)*cos(ky*y)*cos(kz*z) * 300.d0
             u2t = -cos(kx*x)*sin(ky*y)*cos(kz*z) * 300.d0
             u3t = 0.d0

             CALL CKXTY (Xt, IWRK, RWRK, Yt)
             CALL CKRHOY(Pt,Tt,Yt,IWRK,RWRK,rhot)
             call CKUBMS(Tt,Yt,IWRK,RWRK,et)

             cons(i,j,k,irho) = rhot
             cons(i,j,k,imx)  = rhot*u1t
             cons(i,j,k,imy)  = rhot*u2t
             cons(i,j,k,imz)  = rhot*u3t
             cons(i,j,k,iene) = rhot*(et + 0.5d0*(u1t**2 + u2t**2 + u3t**2))

             do n=1,nspecies
                cons(i,j,k,iry1-1+n) = Yt(n)*rhot
             end do

          enddo
       enddo
    enddo
    !$omp end parallel do

  end subroutine init_data_3d

end module init_data_module

max_stars_repo_name: yshoraka/miniSMC

max_issues_repo_name: kiranhegde/miniSMC

max_forks_repo_name: kiranhegde/miniSMC


text:
10 INTERFACE
END INTERFACE
INTERFACE binky
END INTERFACE
END PROGRAM

max_stars_repo_name: ouankou/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab