text:
        SUBROUTINE UNPACK ( *, FILE, A )
      INCLUDE 'PAKBLK.COM'
      INCLUDE 'DSIOF.COM'
      INCLUDE 'XNSTRN.COM'
        COMMON / UNPAKX / ITYPOT, IROBGN, LASROW, INCR
        INTEGER           FILE
        REAL              A(4), IBASE
        DATA              LARGE / 65536 /
        NAME = FILE
        NUM    = NWRDEL( IABS( ITYPOT ) )
        CALL DSIPK1( IBLKD, ITYPOT )
        IF ( IRETRN .EQ. 1 ) GO TO 7700
        IF ( IROBGN .LE. 0 .OR. LASROW .LE. 0 ) GO TO 10
        IROW   = IROBGN
        ILSROW = LASROW
        GO TO 20
10      IROBGN = IBLKD( 4 )
        IROW   = IROBGN
        ILSROW = LARGE
20      CONTINUE
        INDEX2 = 1
        ITYPE  = IBLKD( 13 )
        INDEX1 = ( IBLKD( 5 ) - 1 ) * IBLKD( 14 ) + 1
        NUMINC = NUM * INCR
90      IF ( IBLKD( 4 ) .GT. ILSROW ) GO TO 200
        IF ( ( IBLKD( 4 ) + IBLKD( 6 ) - 1 ) .LT. IROBGN ) GO TO 145
100     IDIFF = ( IBLKD( 4 ) + IBLKD( 7 ) ) - IROW
        IBLKD( 7 ) = IBLKD( 7 ) + 1
        IF ( IDIFF .EQ. 0 ) GO TO 140
        IF ( IDIFF .LT. 0 ) GO TO 142
        DO 130 K = 1, NUM
        DO 110   KKK = 1, IDIFF
        A( INDEX2 + K - 1 + (KKK-1)*NUMINC ) = 0.
110     CONTINUE
130     CONTINUE
        INDEX2 = INDEX2 + IDIFF*NUMINC
        IROW   = IROW + IDIFF
140     IF ( IBLKD(2) .NE. ITYPE ) GO TO 1400
CDIR$ NOVECTOR
        DO 141 K = 1, NUM
        A( INDEX2 + K - 1 ) = IBASE( INDEX1 + K - 1 )
141     CONTINUE
CDIR$ VECTOR
        GO TO 1401
1400    CALL DSUPKC( IBLKD( 2 ), ITYPE, IBASE( INDEX1 ), A( INDEX2 ) )
1401    CONTINUE
        IF ( IROW .GE. ILSROW ) GO TO 225
        IROW = IROW + 1
        INDEX2 = INDEX2 + NUMINC
142     INDEX1 = INDEX1 + IBLKD( 11 )
        IF ( IBLKD( 7 ) .NE. IBLKD( 6 ) ) GO TO 100
145     CALL ENDGET( IBLKD )
        CALL GETSTR( *200, IBLKD )
        INDEX1 = ( IBLKD( 5 ) - 1 ) * IBLKD( 14 ) + 1
        IBLKD( 7 ) = 0
150     IF ( IBLKD( 8 ) .LT. 1 ) GO TO 90
200     IF ( ILSROW .EQ. LARGE ) GO TO 230
        NUMLEF = LASROW - IROW + 1
        IF ( NUMLEF .LE. 0 ) GO TO 225
        DO 220 KK = 1, NUM
        DO 210 K  = 1, NUMLEF
        A( INDEX2 + KK - 1 + (K-1)*NUMINC ) = 0.
210     CONTINUE
220     CONTINUE
225     IF ( IBLKD( 8 ) .GE. 1 ) GO TO 240
        CALL DSSKRC
        CALL DSSDCB
        GO TO 240
230     LASROW = IROW - 1
240     RETURN
7700    RETURN 1
        END


text:
!
!     (c) 2019 Guide Star Engineering, LLC
!     This Software was developed for the US Nuclear Regulatory Commission (US NRC)
!     under contract "Multi-Dimensional Physics Implementation into Fuel Analysis under
!     Steady-state and Transients (FAST)", contract # NRC-HQ-60-17-C-0007
!
!
!    NEMO - Numerical Engine (for) Multiphysics Operators
! Copyright (c) 2007, Stefano Toninel
!                     Gian Marco Bianchi  University of Bologna
!              David P. Schmidt    University of Massachusetts - Amherst
!              Salvatore Filippone University of Rome Tor Vergata
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification,
! are permitted provided that the following conditions are met:
!
!     1. Redistributions of source code must retain the above copyright notice,
!        this list of conditions and the following disclaimer.
!     2. Redistributions in binary form must reproduce the above copyright notice,
!        this list of conditions and the following disclaimer in the documentation
!        and/or other materials provided with the distribution.
!     3. Neither the name of the NEMO project nor the names of its contributors
!        may be used to endorse or promote products derived from this software
!        without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
! ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
! WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
! ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
!---------------------------------------------------------------------------------
!
! $Id: wr_mtx_pattern.f90 3106 2008-05-05 15:07:27Z sfilippo $
!
! Description:
!    Dumps the pattern of a sparse matrix in Matrix Market format, starting
!    from the global C2C adjacency graph.
!
SUBMODULE (tools_output_basics) wr_mtx_pattern_implementation
    IMPLICIT NONE

    CONTAINS

        MODULE PROCEDURE wr_mtx_pattern
            USE class_psblas
            USE class_connectivity
            IMPLICIT NONE
            !
            INTEGER, PARAMETER :: pattern=10
            INTEGER :: i, j, n, ncells, nnz
            INTEGER, POINTER :: iconn(:) => NULL()
            CHARACTER(len=*), PARAMETER :: fmt='(3(i8,1x))'

            CALL sw_out%tic()

            IF(mypnum_() == 0) THEN
                WRITE(*,'(a)',advance='no') ' Now dumping sparsity pattern...'

                OPEN(unit=pattern,file=name)

                ncells = c2c%nel_()

                ! REMARK! The adjacency graph C2C contains only the non-diagonal elements.
                ! One has to add also the elements of the main diagonal in order to get
                ! the right number of non-zero coefficients.
                nnz = c2c%nconn_() + ncells

                WRITE(pattern,fmt) ncells, ncells, nnz

                DO i = 1, ncells

                    CALL c2c%get_ith_conn(iconn,i)
                    n = SIZE(iconn)

                    ! Off-diagonal elements -> 1
                    DO j = 1, n
                        WRITE(pattern,fmt) i, iconn(j), 1
                    END DO

                    ! Main diagonal elements -> -1
                    WRITE(pattern,fmt) i, i, -1
                END DO

                NULLIFY(iconn)

                WRITE(*,*) 'done.'
                WRITE(*,*)

                CLOSE(pattern)
            END IF

            CALL sw_out%toc()

        END PROCEDURE wr_mtx_pattern

END SUBMODULE wr_mtx_pattern_implementation

max_stars_repo_name: sourceryinstitute/MORFEUS-Source

max_issues_repo_name: sourceryinstitute/MORFEUS-Source

max_forks_repo_name: sourceryinstitute/MORFEUS-Source


text:
      program regrid
ccc      parameter(nimax=730,njmax=380,n1max=nimax*2,
      parameter(nimax=2100,njmax=1200,n1max=nimax*2,
     $     niopt=4,nvopt=8)

      real vals(20),ovals(20),args(8)

      real*4 xin(n1max),yin(n1max),xout(n1max),yout(n1max)
      real*4 gxout(n1max),gyout(n1max)
      real*4 glat(n1max),glatb(n1max)

      real*4 dumx(n1max)

      dimension incard(80)
      character*10 carg(niopt)
      character*2 copt(nvopt)

      dimension fld_in(nimax,njmax),area_in(nimax,njmax),
     $     fld_out(nimax,njmax),dum1(nimax,njmax),
     $     area_out(nimax,njmax)

      character card*80,qtitle*24,regrid_method*40
      character incard*1

      logical cyclicxi,wrapxi,gauss_out,
     $     diag_out,vote,bessel,
     $     npole_point,spole_point

      equivalence (incard(1),card)

      data copt/'ig','ba','bl','vt','bs','p1','gg','ma'/
C
C         gfi - input file to the program with the GrADS data object
C
C         gfo - output file to be read by GrADS
C
      iunit_in=8
      iunit_out=10
      iunit_diag=12

      open (iunit_in,file='udf.regrid2.gfi',status='unknown',
     $     form='unformatted',err=810)

      open (iunit_out,file='udf.regrid2.gfo',status='unknown',
     $     form='unformatted')

      open (iunit_diag,file='udf.regrid2.out',
     $     form='formatted',status='unknown')
C
C         constants
C
      pi=3.141592654
      deg2rad=pi/180.0
      rad2deg=180.0/pi
C
C         defaults
C
      diag_out=.false.
      diag_out=.true.
C
C         user defined beginning lat/lon
C
      rlatbeg_set=-999.0
      rlonbeg_set=-999.0
C
C         pole points
C
      npole_point=.false.
      spole_point=.false.
C
C         initialize flag for using calculated gauss lat weights
C
      jsglat=-888
C
C         the default method is box averaging without voting on
C         a unform grid
C
      iregrid_type=1
      iregrid_method=1
      ilinearo=1
      jlinearo=1
      bessel=.false.
      vote=.false.
      gauss_out=.false.
C
C         minimum fractional area covered by data in box averaging
C         to have a defined point
C
      area_min=0.5
C
C         voting defaults --
C         50% of grid box must be covered regardless of the number
C         of candidates
C
          rmin_vote_max=0.5
          rmin_vote_min=0.5
C
C         the first record contains number of args
C         and other parameters for future implementations
C

      read (iunit_in) vals
      nargs=nint(vals(1))
      write(*,*) 'number of arguments = ',nargs
C
C         stop if only one argument
C
      if(nargs.le.1) then
        write(*,12)
 12     format(' ',/
     $       /' ','regrid requires 2 or more arguements:',
     $       /' ',' ')
        go to 999
      endif
C
C         the second record is the field ("expr" in the user defined
C         function table) and has is input using four reads
C
C         grid record 1
C
      read(iunit_in,err=820) vals
C
C         grid parameters
C
      undef=vals(1)
      itype=nint(vals(2))
      jtype=nint(vals(3))
      nii=nint(vals(4))
      nji=nint(vals(5))
      ilineari=nint(vals(6))
      jlineari=nint(vals(7))
      xbeg=vals(8)
C
C         check the input grid type
C
C         abort if not 2-d x,y
C
      if(itype.ne.0.or.jtype.ne.1) go to 830
C
C         get grid increments if uniform
C
      if(ilineari.eq.1) then
        dx=vals(9)
      else
        dx=-999
      endif
      ybeg=vals(10)

      if(jlineari.eq.1) then
        dy=vals(11)
      else
        dy=-999
      endif
C
C         only get time values if a one of the dimensions is time
C
      if(itype.eq.3.or.jtype.eq.3) then
        yrbeg=vals(12)
        mobeg=vals(13)
        dabeg=vals(14)
        hrbeg=vals(15)
        mnbeg=vals(16)
        dmn=vals(17)
        dmo=vals(18)
      endif
C
      write(iunit_diag,102) undef,itype,jtype,nii,nji,
     $     ilineari,jlineari,xbeg,dx,ybeg,dy
 102  format(' ',
     $     /,10x,'          undefined value = ',e13.5,
     $     /,10x,'         x dimension type = ',i3,
     $     /,10x,'         y dimension type = ',i3,
     $     /,10x,'         x dimension size = ',i5,
     $     /,10x,'         y dimension size = ',i5,
     $     /,10x,' uniform grid spacing in x ? ',i2,
     $     /,10x,' uniform grid spacing in y ? ',i2,
     $     /,10x,'   x world coor beginning  = ',e12.4,
     $     /,10x,'                        dx = ',e12.4,
     $     /,10x,'   y world coor beginning  = ',e12.4,
     $     /,10x,'                        dy = ',e12.4)

      if(itype.eq.3.or.jtype.eq.3) then
        write(iunit_diag,104)
     $       yrbeg,mobeg,dabeg,hrbeg,mnbeg,dmn,dmo
 104    format(' ',
     $       /,10x,'                year start = ',e12.4,
     $       /,10x,'               month start = ',e12.4,
     $       /,10x,'                 day start = ',e12.4,
     $       /,10x,'                hour start = ',e12.4,
     $       /,10x,'              minute start = ',e12.4,
     $       /,10x,'           minute increment = ',e12.4,
     $       /,10x,'            month increment = ',e12.4)
      endif

C
C         grid record #2 the data
C
      call read_grid(iunit_in,fld_in,nii,nji,istat)
      if(istat.ne.1) go to 820

      qtitle='input field from GrADS '
      call qprntu(fld_in,qtitle,undef,1,1,nii,nji,4,iunit_diag)
C
C         grid record #3 and 4 -- grid-to-world coordinate maps
C         of the input grid
C
      read(iunit_in,err=820) (xin(i),i=1,nii)
      write(iunit_diag,*) (xin(i),i=1,nii)
      read(iunit_in,err=820) (yin(j),j=1,nji)
      write(iunit_diag,*) (yin(j),j=1,nji)

C
C------------------------------------------------------------
C
C         read arguments after the grid -- two floats
C
C         #1 - dx or # gaussian lats
C         #2 - dy or # gaussian lons
C
C         if only 1 then dx=dy
C         if only 2 then method = box averaging
C
C------------------------------------------------------------
C
C         decrement ngarg by one for the grid
C
      nargs=nargs-1
      if(nargs.ge.2) then
        nargrd=2
      else
        nargrd=1
      endif

      do i=1,nargrd
        read(iunit_in,err=832) args(i)
      end do
C
C         if one argument, then assume uniform lat lon with equal dx
C
      if(nargrd.eq.1) then
        iout_grid_type=1
        iregrid_method=1
        dxout=args(1)
        dyout=args(1)

      else
C
C         if two arguments, dx & dy are specified the
C         we run box ave by default
C
        iout_grid_type=1
        iregrid_method=1
        dxout=args(1)
        dyout=args(2)
C
C         method
C
        iregrid_type=1
        iregrid_method=1
      endif
C
C------------------------------------------------------------
C
C         read the third argument -- the char option string
C
C         this specifies:
C
C         1)  we are telling regrid that the input grid
C         is gaussian and the resolution (number of gaussian latitudes)
C
C         2)  interp method
C
C         3)  starting lon,lat of the (1,1) of the output grid
C
C------------------------------------------------------------
C
C         decrement nargs based on the last reads
C

      nargs=nargs-nargrd
C
C         number of char options
C
      nca=0
C
C         we have char args...
C
      if(nargs.ge.1) then

        read(iunit_in,err=832) incard
        lenic=ichar_len(incard,80)
C
C         parse it and check for proper lengths
C
        nca=1
        ib=1
        do ii=1,lenic

          carg(nca)='          '
          if(incard(ii).eq.'_') then
            ie=ii-1
            ilen=ie-ib+1
            carg(nca)(1:ilen)=card(ib:ie)
            if(carg(nca)(1:2).eq.'ig'.or.carg(nca)(1:2).eq.'ma') then
              if(ilen.eq.2) go to 812
            else
              if(ilen.gt.2) go to 812
            endif

            ib=ie+2
            nca=nca+1

          endif

          if(ii.eq.lenic) then
            ie=lenic
            ilen=ie-ib+1
            carg(nca)(1:ilen)=card(ib:ie)
            if(carg(nca)(1:2).eq.'ig'.or.carg(nca)(1:2).eq.'ma') then
              if(ilen.eq.2) go to 812
            else
              if(ilen.gt.2) go to 812
            endif

          endif

        enddo
C
C         too many args
C
        if(nca.gt.niopt) go to 814
C
C***      process the args
C
        do ii=1,nca
          iok=0
ccc          print*,'ii = ',ii,carg(ii)
          do jj=1,nvopt
            if(carg(ii)(1:2).eq.copt(jj)) then
              iok=1
            endif
          enddo
C
C***      invalid carg
C
          if(iok.ne.1) then
            nca=ii
            go to 812
          endif
C
C***      min area for defined data
C
          if(carg(ii)(1:2).eq.'ma') then
            read(carg(ii)(3:5),'(i3)') i1
            area_min=i1*0.01
C
C         decrease area_min slightly for 100%
C         because of numerical differences between
C         sum of intersections and the actual sfc area
C         of the output grid boxes
C
            if(area_min.ge.1.0) area_min=0.96
          endif
C
C***      number of gaussian lats in input grid
C***      to invoke the more precise calculation
C***      of the box boundaries
C
          if(carg(ii)(1:2).eq.'ig') then

            eps_glat=0.01
            read(carg(ii)(3:5),'(i3)') njglat
C
C         make sure this IS a valid gaussian grid
C         and can be calculated by the routine
C
            if(mod(njglat,4).ne.0.or.njglat.lt.8) then
              write(*,'(/,a,/,a,/)')
     $             'number of input gaussian lats invalid',
     $             '(not a factor of 4) bypass precise weighting'
            else
C
C         get the latitudes and boundaries using the pcmdi routine
C
              call gauss_lat_pcmdi(glat,glatb,njglat)

              do j=1,njglat
                write(*,'(i3,3(2x,f7.3))')
     $               j,glat(j),glatb(j),glatb(j+1)
              end do

              do j=1,nji
                write(*,'(i3,2x,f7.3)') j,yin(j)
              end do
C
C         now see if we can match the input grid to the assumed
C         gaussian grid
C
              jstest=1
              if(yin(jstest).lt.-90.0) jstest=2

              do j=1,njglat
                if(abs(yin(jstest)-glat(j)).le.eps_glat) then
                  if(diag_out) then
                    write(iunit_diag,*)
     $                   ' We have a near match use the calcuated grid'
                    write(iunit_diag,*)
     $                   ' jsglat = ',j,' jstest = ',jstest
                    write(iunit_diag,*)
     $                   ' abs(yin(jstest) =  ',yin(jstest)
                    write(iunit_diag,*)
     $                   ' glat(j) = ',glat(j)
                  endif
                  jsglat=j
C
C         quick exit
C
                  go to 9001
                endif
              end do

 9001         continue

            endif

          endif

        enddo

      endif

C
C
C------------------------------------------------------------
C
C         input grid parameters
C         check for cyclic continuity in x
C
C------------------------------------------------------------
C
      cyclicxi=.false.
      wrapxi=.false.
C
C         make sure the x dimension is longitude and is uniform
C
      if(itype.eq.0.and.ilineari.eq.1) then
        xlen=(nii-1)*dx
        if(xlen.ge.360.0) then
          wrapxi=.true.
          cyclicxi=.true.
          niifix=360.0/dx
          xlen=360.0
        endif

        xlen=nii*dx
        if(xlen.eq.360.0) then
          wrapxi=.false.
          cyclicxi=.true.
        endif

      endif

      if(diag_out)
     $     write(iunit_diag,*) 'cyclicxi: ',cyclicxi,wrapxi
C
C------------------------------------------------------------
C
C         if wrapped in x, then trim the grid
C
C         dump the input data into a dummy array
C         and then load the trimmed grid into the original field array
C
C------------------------------------------------------------
C
      if(wrapxi) then
        call vload2(dum1,fld_in,nii,nji)
        call trim_grid(dum1,nii,niifix,nji,fld_in)
C
C         set the input size to the new trimmed x dimension
C
        nii=niifix
      endif
C
C------------------------------------------------------------
C
C         set up the input grid
C
C------------------------------------------------------------
C
C         boundaries of input grid boxes
C
      rlonbegi=xin(1)
      rlonendi=xin(nii)
      if(cyclicxi) rlonedi=xin(1)+360.0

      rlatbegi=yin(1)
      rlatendi=yin(nji)

      niip1=nii+1
      njip1=nji+1
C
C         x input grid box boundaries
C
      do i=2,nii
        dumx(i)=0.5*(xin(i-1)+xin(i))
      end do

      dumx(1)=xin(1)-0.5*(xin(2)-xin(1))
      dumx(niip1)=xin(nii)+0.5*(xin(nii)-xin(nii-1))

      do i=1,niip1
        xin(i)=dumx(i)
      end do
C
C         y input grid box boundaries
C
      if(jsglat.ne.-888) then
C
C         use the calculated gaussian boundaries
C
        dumx(1)=yin(1)-0.5*(yin(2)-yin(1))
        dumx(njip1)=yin(nji)+0.5*(yin(nji)-yin(nji-1))

        jend=njip1
        if(jend.gt.njglat+1) jend=njglat+1

        if(jstest.eq.1) then

          do j=1,jend
            jj=jsglat+(j-1)
            if(j.ne.1.and.j.ne.njip1)
     $           dumx(j)=0.5*(yin(j-1)+yin(j))
            if(diag_out)
     $           write(iunit_diag,*) 'jstest 1 ',j,jj,dumx(j),glatb(jj)
            dumx(j)=glatb(jj)
          end do

        else if(jstest.eq.2) then

          write(iunit_diag,*) 'jstest 2 1 = ',dumx(1)
          do j=2,jend
            jj=jsglat+(j-jstest)
            if(j.ne.1.and.j.ne.njip1)
     $           dumx(j)=0.5*(yin(j-1)+yin(j))
            if(diag_out)
     $           write(iunit_diag,*) 'jstest 2 ',j,jj,dumx(j),glatb(jj)
            dumx(j)=glatb(jj)
          end do

        end if

      else
C
C         use the input lat -> j map for the boundaries
C
        do j=2,nji
          dumx(j)=0.5*(yin(j-1)+yin(j))
        end do

        dumx(1)=yin(1)-0.5*(yin(2)-yin(1))
        dumx(njip1)=yin(nji)+0.5*(yin(nji)-yin(nji-1))

      endif
C
C         make sure the input grid does not extend beyond the poles
C
      do j=1,njip1
        yin(j)=dumx(j)
        if(yin(j).lt.-90.0) yin(j)=-90.0
        if(yin(j).gt.90.0) yin(j)=90.0
      end do

      if(diag_out) then
        write(iunit_diag,*) 'input grid x boundaries'
        write(iunit_diag,*) (xin(i),i=1,niip1)
        write(iunit_diag,*) 'input grid y boundaries'
        write(iunit_diag,*) (yin(j),j=1,njip1)
      endif

C
C------------------------------------------------------------
C
C         read arguments after the char options
C
C         the first is the grid and the second in the char options
C
C------------------------------------------------------------
C
C         decrement nargs if there were char opts
C
      if(nca.gt.0) nargs=nargs-1

      if(nargs.ne.0) then
        do ii=1,nargs
          read(iunit_in,err=832) args(ii+nargrd)
        enddo
      endif

      istrvt=999
      istrp1=999

      if(nca.ne.0) then

        do ii=1,nca

          if(carg(ii)(1:2).eq.'gg') then

C
C         Gaussian Output Grid
C
            niog=nint(args(1))
            njog=nint(args(2))
            nwaves=niog/3-1
            gauss_out=.true.

C
C         calculate the gaussian latitudes and
C         the grid spacing for the gaussian longitudes
C
            call gauss_lat_nmc(dumx,njog)

            ilinearo=1
            jlinearo=0
            dxout=360.0/float(niog)
            dyout=180.0/float(njog)
            dyoutg=dyout
            iout_grid_type=2

          endif
C
C*****         method
C
          if(carg(ii)(1:2).eq.'vt') then
            iregrid_method=1
            vote=.true.
            if(nargs.gt.0) istrvt=ii
          else if(carg(ii)(1:2).eq.'ba') then
            dxout=args(1)
            dyout=args(2)
            iregrid_method=1
C
C         change the type of map if gauss out grid
C
            if(gauss_out) then
              ilinearo=1
              jlinearo=0
            endif
            vote=.false.
          else if(carg(ii)(1:2).eq.'bl') then
            iregrid_method=2
            bessel=.false.
          else if(carg(ii)(1:2).eq.'bs') then
            iregrid_method=2
            bessel=.true.
          end if
C
C         starting lon/lat of the (1,1) point
C
          if(carg(ii)(1:2).eq.'p1') then
            istrp1=ii
          endif

        end do
C
C         SPECIAL OPTIONS
C
C         case 1 - voting options
C         case 2 - user-defined beginning lat/lon
C
        if(nargs.eq.4) then

          if(istrvt.eq.999.or.istrp1.eq.999) go to 838

          if(istrvt.lt.istrp1) then

            rmin_vote_max=args(3)
            rmin_vote_min=args(4)
            rlonbeg_set=args(5)
            rlatbeg_set=args(6)

          else

            rlonbeg_set=args(3)
            rlatbeg_set=args(4)
            rmin_vote_max=args(5)
            rmin_vote_min=args(6)

          endif

        else if(nargs.eq.2) then
          if(istrvt.ne.999.and.istrp1.ne.999) then
C
C         both vt and p1 have been set, then use defaults for
C         vt and set the start lon,loat
C
            rlonbeg_set=args(3)
            rlatbeg_set=args(4)
          else if(istrvt.ne.999) then
            rmin_vote_max=args(3)
            rmin_vote_min=args(4)
          else if(istrp1.ne.999) then
            rlonbeg_set=args(3)
            rlatbeg_set=args(4)
          else
            go to 839
          endif

        else if(nargs.ne.0) then

          go to 839

        endif

      endif
C
C         QC the input
C
C         #1 - setting the (1,1) point for a gaussian output grid NOT!
C
      if(iout_grid_type.eq.2.and.rlonbeg_set.ne.-999.0) then
        rlonbeg_set=-999.0
        rlatbeg_set=-999.0
      endif

      if(diag_out) then
      write(iunit_diag,*) 'iout_grid_type = ',iout_grid_type
      write(iunit_diag,*) 'iregrid_method = ',iregrid_method
      write(iunit_diag,*) '         dxout = ',dxout
      write(iunit_diag,*) '         dyout = ',dyout
      write(iunit_diag,*) '          vote = ',vote
      write(iunit_diag,*) '        bessel = ',bessel
      write(iunit_diag,*) ' rmin_vote_max = ',rmin_vote_max
      write(iunit_diag,*) ' rmin_vote_min = ',rmin_vote_min
      write(iunit_diag,*) '   rlonbeg_set = ',rlonbeg_set
      write(iunit_diag,*) '   rlatbeg_set = ',rlatbeg_set
      endif

C
C------------------------------------------------------------
C
C         set up the output grid based on the arguments
C
C------------------------------------------------------------
C
C         iout_grid_type 1 - uniform lat/lon in both directions
C
      if(iout_grid_type.eq.1) then

        rlonbego=ifix(rlonbegi/dxout)*dxout
        rlonendo=ifix(rlonendi/dxout)*dxout
        if(cyclicxi) rlonendo=rlonbego+(360.0-dxout)
        rlatbego=ifix(rlatbegi/dyout)*dyout
        rlatendo=ifix(rlatendi/dyout)*dyout
C
C         check for a global input gaussian grid
C
        yleni=yin(njip1)-yin(1)
        if(yleni+dyout.ge.180.0) then
          rlatbego=-90.0
          rlatendo=90.0
        endif

        nio=nint((rlonendo-rlonbego)/dxout)+1
        njo=nint((rlatendo-rlatbego)/dyout)+1

C
C         make sure we have at least a 2x2 output grid
C
        if(nio.lt.2) then
          nio=2
          rlonendo=rlonbego+2.0*dxout
        endif
        if(njo.lt.2) then
          njo=2
          rlatendo=rlonbego+2.0*dyout
        endif
C
C         user defined beginning longitude
C
        if(rlonbeg_set.ne.-999.0) then
          rlonbego=rlonbeg_set
C
C         cylic continuity check
C
          if(cyclicxi) then
            nio=nint(360.0/dxout)
          else
            nio=nint((rlonendi-rlonbego)/dxout)+1
          endif
          if(nio.lt.2) nio=2
          rlonendo=rlonbego+(nio-1)*dxout
        endif

        if(rlatbeg_set.ne.-999.0) then
          rlatbego=rlatbeg_set
          njo=nint((rlatendi-rlatbego)/dyout)+1
          if(njo.lt.2) njo=2
          rlatendo=rlatbego+(njo-1)*dyout
C
C         bounds check
C
          if(rlatendo.gt.90.0) then
            rlatendo=rlatendo-dyout
            njo=njo-1
          endif

        endif

        if(diag_out) then
          write(iunit_diag,*) ' '
          write(iunit_diag,*) 'iout_grid_type=1'
          write(iunit_diag,*)
     $         'i= ',rlonbegi,rlonendi,rlatbegi,rlatendi
          write(iunit_diag,*)
     $         'o= ',rlonbego,rlonendo,rlatbego,rlatendo
          write(iunit_diag,*) 'nio,njo = ',nio,njo
        endif

C
C         type 2 - gaussian in latitude
C
      else if(iout_grid_type.eq.2) then

        rlonbego=ifix(rlonbegi/dxout)*dxout
        rlonendo=ifix(rlonendi/dxout)*dxout
        if(cyclicxi) rlonendo=rlonbego+(360.0-dxout)
        nio=nint((rlonendo-rlonbego))/dxout+1
C
C         limits of the input grid (handles gaussian --> gaussian)
C
        xbegi=yin(1)
        xendi=yin(njip1)

        jj=0
        do j=1,njog
          if(dumx(j).ge.xbegi.and.dumx(j).le.xendi) then
            jj=jj+1
            yout(jj)=dumx(j)
          endif
        end do

        njo=jj
        rlatbego=yout(1)
        rlatendo=yout(njo)

        if(diag_out) then
          write(iunit_diag,*) ' '
          write(iunit_diag,*) 'iout_grid_type=2'
          write(iunit_diag,*) 'dxout = ',dxout
          write(iunit_diag,*)
     $         'i= ',rlonbegi,rlonendi,rlatbegi,rlatendi
          write(iunit_diag,*)
     $         'o= ',rlonbego,rlonendo,rlatbego,rlatendo
          write(iunit_diag,*) 'nio, njo = ',nio,njo
        endif

      else

C
C         invalid grid type; abort
C
        go to 834

      endif
C
C         output grid characteristics to GrADS
C

      if(iregrid_method.eq.1.and..not.vote)
     $     regrid_method='box averaging                  '
      if(iregrid_method.eq.1.and.vote)
     $     regrid_method='box averaging with VOTING      '
      if(iregrid_method.eq.2.and..not.bessel)
     $     regrid_method='bilinear interpolation         '
      if(iregrid_method.eq.2.and.bessel)
     $     regrid_method='bessel interpolation           '

      if(iout_grid_type.eq.1) then
        write(*,'(a)')  ' '
        write(*,'(a)')
     $       'the output grid is UNIFORM lat/lon:'
        write(*,'(a,f5.2,a,f5.2,a)')
     $       'dx = ',dxout,' deg and dy = ',dyout,' deg'
        write(*,'(a,i3,a,i3)')
     $       '# points in i(lon) = ',nio,
     $       '  # points j(lat) = ',njo
        write(*,'(a,f7.2,a,f7.2,a,f6.2,a,f6.2)')
     $       'lon extent = ',rlonbego,' to ',rlonendo,
     $       ' lat extend = ',rlatbego,' to ',rlatendo
      endif

      if(iout_grid_type.eq.2) then
        write(*,'(a)')  ' '
        write(*,'(a,i3,a)')
     $       'the output grid is ~ T',nwaves,' GAUSSIAN:'
        write(*,'(a,f7.2,a,f5.2,a)')
     $       'dx = ',dxout,' deg and dy ~ ',dyoutg,' deg'
        write(*,'(a,i3,a,i3)')
     $       '# points in i(lon) = ',nio,'  # points j(lat) = ',njo
        write(*,'(a,f7.2,a,f7.2,a,f6.2,a,f6.2)')
     $       'lon extent = ',rlonbego,' to ',rlonendo,
     $       '   lat extend = ',rlatbego,' to ',rlatendo
      endif
      write(*,'(a,a)')
     $     'regrid method is: ',regrid_method
      if(vote) then
        write(*,'(a,f4.2,a,f4.2)')
     $       'vote parameters:  max fract area = ',rmin_vote_max,
     $       '  min frac area = ',rmin_vote_min
      endif
      write(*,'(a)')  ' '

C
C         boundaries of the output boxes
C
      niop1=nio+1
      njop1=njo+1

C
C         X
C
      do i=1,niop1
        xout(i)=rlonbego+(i-1)*dxout-0.5*dxout
      end do
C
C         y
C
C         make sure the output grid does not extend beyond the poles
C
      if(iout_grid_type.eq.1) then

        do j=1,njop1
          yout(j)=rlatbego+(j-1)*dyout-0.5*dyout
          if(yout(j).lt.-90.0) yout(j)=-90.0
          if(yout(j).gt.90.0) yout(j)=90.0
        end do
C
C         check for pole points
C
        if(rlatendo.eq.90.0) npole_point=.true.
        if(rlatbego.eq.-90.0) spole_point=.true.

      else if(iout_grid_type.eq.2) then

        dumx(1)=yout(1)-0.5*(yout(2)-yout(1))
        dumx(njop1)=yout(njo)+0.5*(yout(njo)-yout(njo-1))
        do j=2,njo
          dumx(j)=0.5*(yout(j-1)+yout(j))
        end do

        do j=1,njop1
          yout(j)=dumx(j)
        end do

      endif

C
C------------------------------------------------------------
C
C         input-output grid box relationship
C
C------------------------------------------------------------
C
C         make sure longitudes of the input/output grids
C         are in positive deg
C
      if(xin(1).lt.0.0) then
        do i=1,niip1
          xin(i)=xin(i)+360.0
        enddo
      endif

      if(xout(1).lt.0.0) then
        do i=1,niop1
          xout(i)=xout(i)+360.0
        enddo
      endif
C
C         calculate the location of the output grid box boundaries
C         w.r.t. the input grid box boundaries
C
      call in_out_boundaries(xin,yin,nii,nji,
     $     xout,yout,nio,njo,cyclicxi,
     $     niip1,njip1,niop1,njop1,
     $     iunit_diag,
     $     gxout,gyout)

      write(iunit_diag,*) 'output grid x boundaries'
      do i=1,niop1
        write(iunit_diag,*) ' i = ',i,' xout = ',xout(i),
     $       ' gxout = ',gxout(i)
      end do

      write(iunit_diag,*) 'output grid y boundaries'
      do j=1,njop1
        write(iunit_diag,*) ' j = ',j,' yout = ',yout(j),
     $       ' gyout = ',gyout(j)
      end do
C
C         calculate sfc area of each grid box of input grid
C
      call sfc_area(fld_in,xin,yin,undef,nii,nji,
     $     area_in,iunit_diag)

      call sfc_area(fld_out,xout,yout,undef,nio,njo,
     $     area_out,iunit_diag)

cc      qtitle='grid box area on unit s'
cc      call qprntn(area_in,qtitle,1,1,nii,nji,4,iunit_diag)
C
C------------------------------------------------------------
C
C         do the regrid
C
C------------------------------------------------------------
C
C         box averaging or "clumping" with a "voting" option
C         where the output grid equals the value of the
C         input grid which accounts for the most area.  voting
C         is used for discontinuos data such as soil type
C

      if(iregrid_method.eq.1) then

        call box_ave(fld_in,area_in,area_out,area_min,
     $       undef,gxout,gyout,
     $       nii,nji,nio,njo,fld_out,iunit_diag,vote,istat,
     $       rmin_vote_max,rmin_vote_min)
C
C         FNOC bilinear/bessel interpolation
C
      else if(iregrid_method.eq.2) then

        call bssl_interp(fld_in,undef,gxout,gyout,
     $       nii,nji,nio,njo,fld_out,iunit_diag,
     $       cyclicxi,spole_point,npole_point,bessel,istat)

      endif
C
C         check for pole points
C
      if(spole_point.or.npole_point) then
        call fix_poles(fld_out,nio,njo,undef,
     $       spole_point,npole_point)
      endif

      qtitle='fld_out                '
      call qprntu(fld_out,qtitle,undef,1,1,nio,njo,4,iunit_diag)
C
C------------------------------------------------------------
C
C         write out return info for GrADS
C
C------------------------------------------------------------
C
      ovals(1) = 0.0
      write(iunit_out) ovals
C
C         modify the appropriate input grid parameters for
C         defining the output grid for GrADS
C

      vals(4)=float(nio)
      vals(5)=float(njo)

      vals(6)=float(ilinearo)
      vals(7)=float(jlinearo)

      vals(8)=rlonbego
      vals(9)=dxout
      if(iout_grid_type.eq.1) then
        vals(10)=rlatbego
        vals(11)=dyout
      else if(iout_grid_type.eq.2) then
        vals(10)=0.0
        vals(11)=0.0
      endif

C

      write(iunit_out) vals

      call write_grid(iunit_out,fld_out,nio,njo)
C
C         define the lat/lon of the output grid points
C
      if(ilinearo.eq.0) then
        do i=1,nio
          xout(i)=rlonbego+(i-1)*dxout
        end do
        write(iunit_out) (xout(i),i=1,nio)
      endif

      if(iout_grid_type.eq.1) then

        do j=1,njo
          dumx(j)=(yout(j)+yout(j+1))*0.5
          yout(j)=rlatbego+(j-1)*dyout
        end do

      else if(iout_grid_type.eq.2) then

        do j=1,njo
          dumx(j)=(yout(j)+yout(j+1))*0.5
        end do
        do j=1,njo
          yout(j)=dumx(j)
        end do

      endif

      write(iunit_out) (yout(j),j=1,njo)
c
      go to 999
C
C         error conditions
C
 810  continue
      write(*,*) ' '
      write(*,*) 'regrid ERROR 810'
      write(*,*) 'unable to open the GrADS input file to regrid'
      write(*,*) ' '
      go to 999

 812  continue
      write(*,*) ' '
      write(*,*) 'regrid ERROR 812'
      write(*,'(a,a,a)')
     $     'invalid character option in arg #4 = (',
     $     carg(nca)(1:2),') try again'
      write(*,*) ' '
      go to 999

 814  continue
      write(*,*) ' '
      write(*,*) 'regrid ERROR 814'
      write(*,'(a,i1,a)')
     $     'too many options in character arg #4 = (',
     $     nca,') try again'
      write(*,*) ' '
      go to 999

 820  continue
      write(*,*) ' '
      write(*,*) 'regrid ERROR 820'
      write(*,*) 'the first arguement to regrid is not a grid'
      write(*,*) ' '
      go to 999

 830  continue
      write(*,*) ' '
      write(*,*) 'regrid ERROR 830'
      write(*,*) 'the grid must be 2-d and x-y (lon-lat)'
      write(*,*) ' '
      go to 999

 832  continue
      write(*,*) ' '
      write(*,*) 'regrid ERROR 832'
      write(*,*) 'arguments must be numbers'
      write(*,*) ' '
      go to 999

 834  continue
      write(*,*) ' '
      write(*,*) 'regrid ERROR 834'
      write(*,*) 'invalid output grid iout_grid_type  = ',
     $     iout_grid_type
      write(*,*) ' '
      go to 999

 836  continue
      write(*,*) ' '
      write(*,*) 'regrid ERROR 836'
      write(*,*) 'invalid interpolation option',
     $     ' iregrid_type = ',iregrid_type
      write(*,*) ' '
      go to 999

 838  continue
      write(*,*) ' '
      write(*,*) 'regrid ERROR 838'
      write(*,*) '4 or more special case arguments, not a special case'
      write(*,*) ' '
      go to 999

 839  continue
      write(*,*) ' '
      write(*,*) 'regrid ERROR 839'
      write(*,*) 'too few or wrong number of options for special cases'
      write(*,*) ' '
      go to 999

 840  continue
      write(*,*) ' '
      write(*,*) 'regrid ERROR 840'
      write(*,*) 'rmin_vote outside the interval [1,0]'
      write(*,*) ' '
      go to 999

 842  continue
      write(*,*) ' '
      write(*,*) 'regrid ERROR 842'
      write(*,*) 'dlon/dlat greater than zero when attempting'
      write(*,*) 'to set beginning lat/lon of the output grid'
      write(*,*) ' '
      go to 999

 999  continue
      stop
      end

      subroutine box_ave(fld_in,area_in,area_out,area_min,
     $     undef,gxout,gyout,
     $     nii,nji,nio,njo,fld_out,iunit_diag,vote,istat,
     $     rmin_vote_max,rmin_vote_min)

      parameter (na=300)
      dimension fld_in(nii,nji),area_in(nii,nji),
     $     fld_out(nio,njo),area_out(nio,njo),
     $     gxout(nio+1),gyout(njo+1)

      dimension area_box(na),fld_box(na),dxdy_box(na),ifld_rank(na),
     $     fld_cand(na),area_cand(na),dxdy_cand(na),rmin_dxdy_vote(3)

      logical cyclicx,diag_out,vote

      diag_out=.false.
ccc      diag_out=.true.
      istat=1
C
C         tolerance for checking whether grid is undefined (zero sfc area)
C
      eps=1.0e-12
C
C         minimum fractional area of the output grid box
C         in order to have a winner in the voting
C
      rmin_dxdy_vote(1)=rmin_vote_max
      rmin_dxdy_vote(2)=(rmin_vote_max+rmin_vote_min)*0.5
      rmin_dxdy_vote(3)=rmin_vote_min

      ibb=1
      iee=nio
      jbb=1
      jee=njo

      do j=jbb,jee
        do i=ibb,iee

          ib=int(gxout(i))+1
          ie=int(gxout(i+1))+1
          jb=int(gyout(j))+1
          je=int(gyout(j+1))+1
C
C         check for exceeding n pole
C
          if(je.gt.nji) je=nji

C
C         cyclic continuity in x
C
          cyclicx=.false.
          if(ie.lt.ib) then
            ie=ie+nii
            cyclicx=.true.
          endif

          if(diag_out) then
            write(iunit_diag,*) ' '
            write(iunit_diag,*)
     $           'i,j,ib,ie,jb,je = ',i,j,ib,ie,jb,je
          endif
C
C         initialize the counter for intersecting grid boxes
C
          icnt=0
C
C         CASE 1:  only one input grid box in output grid box
C
          if(ib.eq.ie.and.jb.eq.je) then

            icnt=1
            dxdy_box(icnt)=1.0
            area_box(icnt)=area_in(ib,jb)
            if(area_box(icnt).eq.0.0) dxdy_box(icnt)=0.0
            fld_box(icnt)=fld_in(ib,jb)

            if(diag_out) then
              write(iunit_diag,*)
     $             'ib=ie and jb=je ',area_box(icnt),dxdy_box(icnt),
     $             fld_box(icnt)
            endif

          else if(ib.eq.ie) then
C
C         CASE 2:  intersecting boxes in y only
C
            ii=ib
            dx=gxout(i+1)-gxout(i)
            dxout=1.0

            do jj=jb,je

              icnt=icnt+1

              if(icnt.gt.na) then
                istat=0
                go to 999
              endif

              if(jj.eq.jb) then
                dy=float(jj)-gyout(j)
              else if(jj.eq.je) then
                dy=gyout(j+1)-float(jj-1)
              else
                dy=1.0
              endif

              if(jj.eq.jb.or.jj.eq.je) then
                dyout=dy/(gyout(j+1)-gyout(j))
              else
                dyout=1.0
              endif

              dxdy_box(icnt)=dxout*dyout
              area_box(icnt)=dx*dy*area_in(ii,jj)
              if(area_in(ii,jj).eq.0.0) dxdy_box(icnt)=0.0
              fld_box(icnt)=fld_in(ii,jj)

              if(diag_out) then
                write(iunit_diag,*)
     $               'ib=ie ',ii,jj,dx,dy,area_in(ii,jj)
                write(iunit_diag,*)
     $               'area_box ... ',area_box(icnt),fld_box(icnt)
                write(iunit_diag,*)
     $               'dx dyout ... ',dxout,dyout,dxdy_box(icnt)
              endif

            end do

          else if(jb.eq.je) then
C
C         CASE 3:  intersecting boxes in x only
C
            jj=jb
            dy=gyout(j+1)-gyout(j)
            dyout=1.0

            do ii=ib,ie

              icnt=icnt+1
              if(icnt.gt.na) then
                istat=0
                go to 999
              endif

              ii0=ii
              if(cyclicx.and.ii0.gt.nii) ii0=ii-nii
              if(ii.eq.ib) then
                dx=float(ii)-gxout(i)
              else if(ii.eq.ie) then
                x0=float(ii-1)
                if(cyclicx) x0=float(ii0)-1.0
                dx=gxout(i+1)-x0
              else
                dx=1.0
              endif

              if(ii.eq.ib.or.ii.eq.ie) then
                dxout=dx/(gxout(i+1)-gxout(i))
              else
                dxout=1.0
              endif

              dxdy_box(icnt)=dxout*dyout

              area_box(icnt)=dx*dy*area_in(ii0,jj)
              if(area_in(ii0,jj).eq.0.0) dxdy_box(icnt)=0.0
              fld_box(icnt)=fld_in(ii0,jj)

              if(diag_out) then
                write(iunit_diag,*)
     $               'jb=je ',ii,ii0,jj,dx,dy,area_in(ii0,jj)
                write(iunit_diag,*)
     $               'area_box ... ',area_box(icnt),fld_box(icnt)
                write(iunit_diag,*)
     $               'dx dyout ... ',dxout,dyout,dxdy_box(icnt)
              endif

            end do

          else
C
C         CASE 4:  intersecting boxes in both directions
C
            do jj=jb,je

              if(jj.eq.jb) then
                dy=float(jj)-gyout(j)
              else if(jj.eq.je) then
                dy=gyout(j+1)-float(jj-1)
              else
                dy=1.0
              endif

              if(jj.eq.jb.or.jj.eq.je) then
                dyout=dy/(gyout(j+1)-gyout(j))
              else
                dyout=1.0
              endif

              do ii=ib,ie
                icnt=icnt+1
                if(icnt.gt.na) then
                  istat=0
                  go to 999
                endif

                ii0=ii
                if(cyclicx.and.ii0.gt.nii) ii0=ii-nii
                if(ii.eq.ib) then
                  dx=float(ii)-gxout(i)
                else if(ii.eq.ie) then
                  x0=float(ii-1)
                  if(cyclicx) x0=float(ii0)-1.0
                  dx=gxout(i+1)-x0
                else
                  dx=1.0
                endif

                if(ii.eq.ib.or.ii.eq.ie) then
                  dxout=dx/(gxout(i+1)-gxout(i))
                else
                  dxout=1.0
                endif

                dxdy_box(icnt)=dxout*dyout
                area_box(icnt)=dx*dy*area_in(ii0,jj)
                if(area_in(ii0,jj).eq.0.0) dxdy_box(icnt)=0.0
                fld_box(icnt)=fld_in(ii0,jj)
                if(diag_out) then
                  write(iunit_diag,*) 'ib.ne.ib.and.jb.ne.je ',
     $                 ii,ii0,jj,dx,dy,area_in(ii0,jj)
                  write(iunit_diag,*)
     $                 'area_box ... ',area_box(icnt),fld_box(icnt)
                  write(iunit_diag,*)
     $                 'dx dyout ... ',dxout,dyout,dxdy_box(icnt)
                endif

              end do

            enddo

          endif
C
C         integrate or vote for the average value
C
          if(vote) then
C
C         voting routine; first get total area
C

            tot_area=0.0
            do ii=1,icnt
              tot_area=tot_area+area_box(ii)
            end do

            if(diag_out) then
              write(iunit_diag,*) 'i,j,icnt,tot_area = ',
     $             i,j,icnt,tot_area
            endif

            if(tot_area.le.eps) then
              fld_out(i,j)=undef
              go to 100
            endif

            if(icnt.eq.1) then
C
C         USSR election -- only one "candidate" to vote for
C
C         check if the the total area is greater
C         than the minimum required to hold the election (e.g., 0.5)
C
              if(diag_out) then
                write(iunit_diag,*)
     $               'USSR:  ',dxdy_box(1),rmin_dxdy_vote(1),fld_box(1)
              endif

              if(dxdy_box(1).lt.rmin_dxdy_vote(1)) then
                fld_out(i,j)=undef
              else
                fld_out(i,j)=fld_box(1)
              endif

              go to 100

            else if(icnt.eq.2) then
C
C         USA election -- two-party, two-candidate race; area wins
C
              if(diag_out) then
                write(iunit_diag,*)'USA: ',dxdy_box(1),dxdy_box(2),
     $               rmin_vote
              endif

              if(dxdy_box(1).eq.0.0.or.dxdy_box(2).eq.0.0) then
                rmin_vote=rmin_dxdy_vote(1)
              else if(fld_box(1).eq.fld_box(2)) then
                rmin_vote=rmin_dxdy_vote(1)
              else
                rmin_vote=rmin_dxdy_vote(2)
              endif

              if((dxdy_box(1).ge.dxdy_box(2)).and.
     $             (dxdy_box(1).gt.rmin_vote)) then
                fld_out(i,j)=fld_box(1)
              else if(dxdy_box(2).gt.rmin_vote) then
                fld_out(i,j)=fld_box(2)
C
C         case where both candidates are the same in the two-person race
C
              else if(fld_box(1).eq.fld_box(2).and.
     $               ((dxdy_box(1)+dxdy_box(2)).gt.rmin_vote)) then
                fld_out(i,j)=fld_box(2)
              else
                fld_out(i,j)=undef
              endif

            else
C
C         a wide open election - three or more candidates
C
C         sort the data by surface area using
C         the numercial recipes routine --  indexx
C
              call indexx(icnt,fld_box,ifld_rank)

              if(diag_out) then
                write(iunit_diag,*)
     $               'fld_box = ',(fld_box(ipp),ipp=1,icnt)
                write(iunit_diag,*)
     $               'area_box = ',(area_box(ipp),ipp=1,icnt)
                write(iunit_diag,*)
     $               'ifld_rank = ',(ifld_rank(ipp),ipp=1,icnt)
                write(iunit_diag,*)
              endif
C
C         the indexes are in reverse order, with the
C         biggest fld element in the last element of the array
C         first check if the biggest is in the majority
C
C         set up the candidates
C
              ncand=1
              it1=ifld_rank(1)

              area_cand(ncand)=area_box(it1)
              fld_cand(ncand)=fld_box(it1)
              dxdy_cand(ncand)=dxdy_box(it1)

              do ii=2,icnt

                i1=ii-1
                i2=ii
                it1=ifld_rank(i1)
                it2=ifld_rank(i2)

                if(fld_box(it1).eq.fld_box(it2)) then
                  area_cand(ncand)=area_cand(ncand)+area_box(it2)
                  dxdy_cand(ncand)=dxdy_cand(ncand)+dxdy_box(it2)
                else
                  ncand=ncand+1
                  area_cand(ncand)=area_box(it2)
                  dxdy_cand(ncand)=dxdy_box(it2)
                  fld_cand(ncand)=fld_box(it2)
                endif

              end do
              if(diag_out) then
                write(iunit_diag,*) 'ncand = ',ncand
                write(iunit_diag,*)
     $               'fld_cand = ',(fld_cand(ipp),ipp=1,ncand)
                write(iunit_diag,*)
     $               'area_cand = ',(area_cand(ipp),ipp=1,ncand)
                write(iunit_diag,*)
     $               'dxdy_cand = ',(dxdy_cand(ipp),ipp=1,ncand)
              endif
C
C         if one candidate, all done
C
              if(ncand.eq.1) then
                if(dxdy_cand(1).gt.rmin_dxdy_vote(3)) then
                  fld_out(i,j)=fld_cand(1)
                else
                  fld_out(i,j)=undef
                endif
                go to 100
              else

C
C         the candidate with the most area is the winner
C
                area_max=0.0
                do ii=1,ncand
                  if(area_cand(ii).gt.area_max) then
                    iamax=ii
                    area_max=area_cand(ii)
                  endif
                end do

                if(ncand.le.2) then
                  rmin_vote=rmin_dxdy_vote(ncand)
                else
                  rmin_vote=rmin_dxdy_vote(3)
                endif

                if(dxdy_cand(iamax).gt.rmin_vote) then
                  fld_out(i,j)=fld_cand(iamax)
                else
                  fld_out(i,j)=undef
                endif

                if(diag_out) then
                  write(iunit_diag,*) ' '
                  write(iunit_diag,*) 'the winner is...',
     $                 iamax,area_max,dxdy_cand(iamax),
     $                 rmin_vote,fld_out(i,j)
                endif

                go to 100

              endif

            endif

          else
C
C         area integrate
C
            tot_fld=0.0
            tot_area=0.0
C
            do ii=1,icnt
              tot_fld=tot_fld+fld_box(ii)*area_box(ii)
              tot_area=tot_area+area_box(ii)
            end do

            if(tot_area.gt.area_out(i,j)*area_min) then
              fld_out(i,j)=tot_fld/tot_area
            else
              fld_out(i,j)=undef
            endif

            if(diag_out) then
              write(iunit_diag,*) 'qqq ',tot_area,' ',
     $             area_out(i,j)*area_min,' ',
     $             area_out(i,j),' ',fld_out(i,j)
            endif

          endif

 100      continue

          if(diag_out) then
            write(iunit_diag,*)
     $           'i,j fld_out(i,j) = ',i,j,fld_out(i,j)
          endif

        end do
      end do

 999  continue

      return
      end

      subroutine bssl_interp(fld_in,undef,gxout,gyout,
     $     nii,nji,nio,njo,fld_out,iunit_diag,
     $     cyclicxi,spole_point,npole_point,bessel,istat)

      dimension fld_in(nii,nji),fld_out(nio,njo),
     $     gxout(nio+1),gyout(njo+1)

      dimension fr(4)

      logical cyclicxi,diag_out,bessel,spole_point,npole_point

      diag_out=.false.
ccc      diag_out=.true.
      jchk=2
      istat=1
C
C         convert the box boundaries to grid point center
C
      do i=1,nio
C
C         check for crossing the boundaries
C         the only way for this to occur is if the field is
C         cyclically continuous in x
C

        if(gxout(i+1).lt.gxout(i)) then
          gxout(i)=((gxout(i)-float(nii))+gxout(i+1))*0.5+0.5
          if(gxout(i).lt.1.0) gxout(i)=float(nii)+gxout(i)
        else
          gxout(i)=(gxout(i)+gxout(i+1))*0.5+0.5
          if(gxout(i).lt.1.0) gxout(i)=float(nii)+gxout(i)
        endif
      end do

      do j=1,njo

        gyout(j)=(gyout(j)+gyout(j+1))*0.5+0.5
C
C         check if a pole points on the input grid
C
        if(spole_point.and.gyout(j).lt.1.0) gyout(j)=1.0
        if(npole_point.and.gyout(j).gt.nji+0.5) gyout(j)=nji

      end do

      ibb=1
      iee=nio
      jbb=1
      jee=njo

      niim1=nii-1
      njim1=nji-1

      do j=jbb,jee
        do i=ibb,iee

          ic=ifix(gxout(i))
          jc=ifix(gyout(j))

          icp1=ic+1
          if(cyclicxi.and.icp1.gt.nii) icp1=icp1-nii
          jcp1=jc+1

          if(diag_out.and.j.eq.jchk) write(iunit_diag,*)
     $         'i,j,gxout,gyout,ic,jc',
     $         i,j,gxout(i),gyout(j),ic,jc

          if((jc.lt.1.or.jc.gt.nji).or.
     $         (.not.cyclicxi.and.(ic.lt.1.or.ic.gt.nii))) then
            fld_out(i,j)=undef
            if(diag_out.and.j.eq.jchk)
     $           write(iunit_diag,*) 'out of bounds '
            go to 100
          end if
C
C------------------------------------------------
C
C         bilinear/bessel interpolation based on the FNOC routine
C         bssl5 by D. Hensen, FNOC
C
C------------------------------------------------
C
          r=gxout(i)-ic
          s=gyout(j)-jc
C
C         interior check
C
          if((jc.ge.2.and.jc.lt.njim1.and.cyclicxi).or.
     $         (ic.ge.2.and.jc.ge.2.and.
     $         ic.lt.niim1.and.jc.lt.njim1)) go to 10
C
C         border zone check
C
          if((jc.lt.nji.and.cyclicxi).or.
     $         (ic.lt.nii.and.jc.lt.nji)) go to 5
C
C------------------------------------------------
C
C         top and right edge processing
C
C------------------------------------------------
C
          if(ic.eq.nii.and.jc.eq.nji) then

            fld_out(i,j)=fld_in(nii,nji)
            if(diag_out.and.j.eq.jchk)
     $           write(iunit_diag,*) 'upper right hand corenr'

          else if(ic.eq.nii) then

            if(diag_out.and.j.eq.jchk)
     $           write(iunit_diag,*) 'right edge'

            if(fld_in(ic,jc).ne.undef.and.
     $           fld_in(ic,jcp1).ne.undef) then
              fld_out(i,j)=(1.0-s)*fld_in(ic,jc)+
     $             s*fld_in(ic,jcp1)
            else
              fld_out(i,j)=undef
            endif

          else if(jc.eq.nji) then

            if(diag_out.and.j.eq.jchk)
     $           write(iunit_diag,*) 'top edge'

            if(fld_in(ic,jc).ne.undef.and.
     $           fld_in(icp1,jc).ne.undef) then
              fld_out(i,j)=(1.0-r)*fld_in(ic,jc)+
     $             r*fld_in(icp1,jc)
            else
              fld_out(i,j)=undef
            endif

          endif

          go to 100

 5        continue
C
C------------------------------------------------
C
C         border zone; bilinear
C
C------------------------------------------------
C
          iok_bilinear=1
          if(fld_in(ic,jc).eq.undef.or.
     $         fld_in(icp1,jc).eq.undef.or.
     $         fld_in(ic,jcp1).eq.undef.or.
     $         fld_in(icp1,jcp1).eq.undef) iok_bilinear=0

          if(diag_out.and.j.eq.jchk)
     $         write(iunit_diag,*)
     $         'border zone, iok_bilinear= ',iok_bilinear

          if(iok_bilinear.eq.0) then
            fld_out(i,j)=undef
          else
            fld_out(i,j)=(1.-s)*((1.-r)*fld_in(ic,jc)
     $           +r*fld_in(icp1,jc))
     $           +s*((1.-r)*fld_in(ic,jcp1)
     $           +r*fld_in(icp1,jcp1))
          endif

          go to 100

 10       continue
C
C------------------------------------------------
C
C         interior zone
C
C------------------------------------------------
C
C         first check if bilinear is OK
C
          iok_bilinear=1
          if(fld_in(ic,jc).eq.undef.or.
     $         fld_in(icp1,jc).eq.undef.or.
     $         fld_in(ic,jcp1).eq.undef.or.
     $         fld_in(icp1,jcp1).eq.undef) iok_bilinear=0

          if(diag_out.and.j.eq.jchk)
     $         write(iunit_diag,*)
     $         'interior zone, iok_bilinear= ',iok_bilinear

          if(iok_bilinear.eq.0) then
            fld_out(i,j)=undef
            go to 100
          else
C
C         bilinear value is the first guess
C
            fld_out(i,j)=(1.-s)*((1.-r)*fld_in(ic,jc)
     $           +r*fld_in(icp1,jc))
     $           +s*((1.-r)*fld_in(ic,jcp1)
     $           +r*fld_in(icp1,jcp1))
C
C         exit if only doing bilinear
C
            if(.not.bessel) go to 100

          endif

C
C         interpolate 4 columns (i-1,i,i+1,i+2) to j+s and store in fr(1)
C         through fr(4)
C
          r1=r-0.5
          r2=r*(r-1.)*0.5
          r3=r1*r2*0.3333333333334
          s1=s-0.5
          s2=s*(s-1.)*0.5
          s3=s1*s2*0.3333333333334
C
          k=0
          im1=ic-1
          ip2=ic+2

          if(diag_out.and.j.eq.jck) write(iunit_diag,*) 'bessel interp'

C
          do ii=im1,ip2

            k=k+1

            i1=ii

            if(cyclicxi.and.i1.lt.1) i1=nii-i1
            if(cyclicxi.and.i1.gt.nii) i1=i1-nii

            j1=jc

            j1p1=j1+1
            j1p2=j1+2
            j1m1=j1-1

            fij=fld_in(i1,j1)
            fijp1=fld_in(i1,j1p1)
            fijp2=fld_in(i1,j1p2)
            fijm1=fld_in(i1,j1m1)
            if(diag_out.and.j.eq.jchk)
     $           write(iunit_diag,*) 'i1,j1 = ',i1,j1,
     $           fij,fijp1,fijp2,fijm1
C
C         exit if any value undefined
C
            if(fij.eq.undef.or.
     $           fijp1.eq.undef.or.
     $           fijp2.eq.undef.or.
     $           fijm1.eq.undef) go to 100

            u=(fij+fijp1)*0.5
            del=fijp1-fij
            udel2=(fijp2-fijp1+fijm1-fij)*0.5
            del3=fijp2-fijp1-2.*del+fij-fijm1

            fr(k)=u+s1*del+s2*udel2+s3*del3

          end do

C
C         interpolate the fr row to ii+r
C
          u=(fr(2)+fr(3))*0.5
          del=fr(3)-fr(2)
          udel2=(fr(4)-fr(3)+fr(1)-fr(2))*0.5
          del3=fr(4)-fr(3)-2.*del+fr(2)-fr(1)

          fld_out(i,j)=u+r1*del+r2*udel2+r3*del3

 100      continue
          if(diag_out.and.j.eq.jchk)
     $         write(iunit_diag,*) 'interp value = ',fld_out(i,j)

        end do
      end do

 999  continue

      return
      end

      subroutine indexx(n,arrin,indx)
      dimension arrin(n),indx(n)
      do 11 j=1,n
        indx(j)=j
11    continue
      l=n/2+1
      ir=n
10    continue
        if(l.gt.1)then
          l=l-1
          indxt=indx(l)
          q=arrin(indxt)
        else
          indxt=indx(ir)
          q=arrin(indxt)
          indx(ir)=indx(1)
          ir=ir-1
          if(ir.eq.1)then
            indx(1)=indxt
            return
          endif
        endif
        i=l
        j=l+l
20      if(j.le.ir)then
          if(j.lt.ir)then
            if(arrin(indx(j)).lt.arrin(indx(j+1)))j=j+1
          endif
          if(q.lt.arrin(indx(j)))then
            indx(i)=indx(j)
            i=j
            j=j+j
          else
            j=ir+1
          endif
        go to 20
        endif
        indx(i)=indxt
      go to 10
      end

      subroutine in_out_boundaries(xin,yin,nii,nji,
     $     xout,yout,nio,njo,cyclicxi,
     $     niip1,njip1,niop1,njop1,
     $     iunit_diag,
     $     gxout,gyout)
C
C  Purpose:
C
C         calculate the location of grid box boundaries of
C         an "output" grid w.r.t. an "input" grid
C
C         used in a 2-D regriding process, i.e.,
C         input --> output
C
C  Input variables:
C
C         xin - longitudes of the input grid
C         yin - latitudes of the input grid
C         xout - longtitudes of the output grid
C         yout - latitudes of the output grid
C
C         cyclicxi - flag whether the input grid is cyclically continuous
C         in x
C
C         nii - size of the x dimension of the input grid
C         nji - size of the y dimension of the input grid
C         nio - size of the x dimension of the output grid
C         njo - size of the y dimension of the output grid
C
C         niip1= nii+1, etc.
C
C         iunit_diag - unit number to write diagnositics
C
C  Output variables:
C
C         gxout - x location of the output grid in input grid units
C         gyout - y location of the output grid in input grid units
C

      dimension xin(*),yin(*),xout(*),yout(*),
     $     gxout(*),gyout(*)

      logical cyclicxi

C
C         locate the output grid w.r.t. input grid in x
C         ALLOW FOR CYCLIC CONTINUITY IN X!!
C
      ibeg=1
      iend=niip1

      do i=1,niop1

        x0=xout(i)

        do ii=ibeg,iend
C
C         point is before start of input grid
C
          x1=xin(ii)
          x2=xin(niop1)

          if(x0.lt.x1) then
            if(cyclicxi) then
              ic=nii
              do while(x0.lt.x1)
                x2=x1
                x1=xin(ic)-360.0
                ic=ic-1
              enddo
              dx0=(x0-x1)/(x2-x1)
              gxout(i)=float(ic)+dx0
            else
              gxout(i)=0.0
            endif

            ibeg=1
            go to 50

C
C         point is within the input grid
C
          else if(ii.ge.1.and.ii.le.nii) then

            x2=xin(ii+1)

cccc            write(iunit_diag,*) ibeg,i,ii,x0,x1,x2

            if(x0.ge.x1.and.x0.lt.x2) then
              dx0=(x0-x1)/(x2-x1)
              gxout(i)=float(ii-1)+dx0
              ibeg=ii
              go to 50
            endif

C
C         point is beyond input grid
C
          else if(x0.ge.x2) then
            ic=2
            if(cyclicxi) then
              x1=xin(niip1)
              x2=xin(ic)+360.0

cccc              write(iunit_diag,*) 'after',i,ic,x0,x1,x2
              do while(x0.lt.x1.or.x0.ge.x2)
                ic=ic+1
                x1=x2
                x2=xin(ic)+360.0
              end do
              dx0=(x0-x1)/(x2-x1)
              gxout(i)=float(ic-2)+dx0
            else
              gxout(i)=float(nii)
            endif

            ibeg=niip1
            go to 50

          endif

        end do

 50     continue

      end do
C
C         locate the output grid w.r.t. input grid in y
C         NO CYCLIC CONTINUITY!!
C
      jbeg=1
      jend=njip1

      do j=1,njop1

        y0=yout(j)

        do jj=jbeg,jend
C
C         point is before start of input grid
C
          y1=yin(jj)

          if(y0.lt.y1) then

            gyout(j)=0.0
            jbeg=1
cccc            write(iunit_diag,*) 'before y ',j,jj,gyout(j)
            go to 60
C
C         point is within the input grid
C
          else if(jj.ge.1.and.jj.le.nji) then

            y2=yin(jj+1)

cccc            write(iunit_diag,*) 'jbeg... ',jbeg,j,jj,y0,y1,y2

            if(y0.ge.y1.and.y0.lt.y2) then
              dy0=(y0-y1)/(y2-y1)
              gyout(j)=float(jj-1)+dy0
              jbeg=jj
              go to 60
            endif
C
C         point is beyond input grid
C
          else if(y0.ge.y2) then
            gyout(j)=float(nji)
            jbeg=njip1
cccc            write(iunit_diag,*) 'after y ',j,jj,gyout(j)
            go to 60
          endif

        end do

 60     continue

      end do

      return
      end

      subroutine avg_pole(a,m,n)
C
Csss      routine to replace pole value with average
Csss      at the penultimate point
C
      dimension a(m,n)
      ave_2=0.0
      ave_nm1=0.0
      do i=1,m
        ave_2=ave_2+a(i,2)
        ave_nm1=ave_nm1+a(i,n-1)
      end do
      ave_2=ave_2/m
      ave_nm1=ave_nm1/m
      do i=1,m
        a(i,1)=ave_2
        a(i,n)=ave_nm1
      end do
      return
      end

      subroutine bsslz1(bes,n)
c
c
      implicit real*8(a-h,o-z)
      dimension bes(n)
      dimension bz(50)
c
      data pi/3.14159265358979d0/
      data bz  / 2.4048255577d0, 5.5200781103d0,
     $  8.6537279129d0,11.7915344391d0,14.9309177086d0,18.0710639679d0,
     $ 21.2116366299d0,24.3524715308d0,27.4934791320d0,30.6346064684d0,
     $ 33.7758202136d0,36.9170983537d0,40.0584257646d0,43.1997917132d0,
     $ 46.3411883717d0,49.4826098974d0,52.6240518411d0,55.7655107550d0,
     $ 58.9069839261d0,62.0484691902d0,65.1899648002d0,68.3314693299d0,
     $ 71.4729816036d0,74.6145006437d0,77.7560256304d0,80.8975558711d0,
     $ 84.0390907769d0,87.1806298436d0,90.3221726372d0,93.4637187819d0,
     $ 96.6052679510d0,99.7468198587d0,102.888374254d0,106.029930916d0,
     $ 109.171489649d0,112.313050280d0,115.454612653d0,118.596176630d0,
     $ 121.737742088d0,124.879308913d0,128.020877005d0,131.162446275d0,
     $ 134.304016638d0,137.445588020d0,140.587160352d0,143.728733573d0,
     $ 146.870307625d0,150.011882457d0,153.153458019d0,156.295034268d0/

      nn=n
      if(n.le.50) go to 12
      bes(50)=bz(50)
      do 5 j=51,n
    5 bes(j)=bes(j-1)+pi
      nn=49
   12 do 15 j=1,nn
   15 bes(j)=bz(j)
      return
      end
      function esatw(t)
c
c     t is temperature of air in deg celcius.
c     ta is temperature in deg kelvin
c     p is pressure in mb
c     esatw is saturation vapor pressure in mb (over water)
c
      data ps/1013.246/,ts/373.16/
c
      ta = t
      if(t .lt. 100.0) then
         ta = t + 273.16
      end if
      e1=11.344*(1.0-ta/ts)
      e2=-3.49149*(ts/ta-1.0)
      f1=-7.90298*(ts/ta-1.0)
      f2=5.02808*alog10(ts/ta)
      f3=-1.3816*(10.0**e1-1.0)*1.e-7
      f4=8.1328*(10.0**e2-1.0)*1.e-3
      f5=alog10(ps)
      f=f1+f2+f3+f4+f5
      esatw=10.0**f
      return
      end
c
      function esati(t)
c
c     t is temperature of air in deg celcius.
c     ta is temperature in deg kelvin
c     p is pressure in mb
c     esati is saturation vapor pressure with respect to ice in mb
c
      data p0/6.1071/,t0/273.16/
c
      t = ta
      if(t .lt. 100.0) then
        ta = t + 273.16
      end if
      f1=-9.09718*(t0/ta-1.0)
      f2=-3.56654*alog10(t0/ta)
      f3=0.876793*(1.0-ta/t0)
      f4=alog10(p0)
      esati = 10.0**(f1+f2+f3+f4)
      return
      end

      subroutine fix_grid(a,ni,nj)
      dimension a(ni,nj)
      do i=1,ni
        do j=1,nj
          a(i,j)=i+j
        end do
      end do
      return
      end

      subroutine fix_poles(fld_out,nio,njo,undef,
     $     spole_point,npole_point)

      dimension fld_out(nio,njo)
      logical spole_point,npole_point

      rmeans=0.0
      rmeann=0.0
      icnts=0
      icntn=0

      do i=1,nio

        if(fld_out(i,1).ne.undef) then
          rmeans=rmeans+fld_out(i,1)
          icnts=icnts+1
        endif

        if(fld_out(i,njo).ne.undef) then
          rmeann=rmeann+fld_out(i,njo)
          icntn=icntn+1
        endif

      end do

      if(icnts.gt.0) then
        rmeans=rmeans/float(icnts)
      else
        rmeans=undef
      endif

      if(icntn.gt.0) then
        rmeann=rmeann/float(icntn)
      else
        rmeann=undef
      endif

      if(spole_point) then
        do i=1,nio
          fld_out(i,1)=rmeans
        end do
      end if

      if(npole_point) then
        do i=1,nio
          fld_out(i,njo)=rmeann
        end do
      end if

      return
      end

c
c$$$  subprogram documentation block
c                .      .    .                                       .
c subprogram:    gaulat      calculates gaussian grid latitudes
c   prgmmr: s. j. lord       org: w/nmc22    date: 91-06-06
c
c abstract: calculates gaussian grid latitudes
c
c program history log:
c   91-06-06  s. j. lord - copied from kanamitsu library
c   930921 m.fiorino - changed from colatitude to latitude
c   yy-mm-dd  modifier1   description of change
c   yy-mm-dd  modifier2   description of change
c
c usage:    call pgm-name(inarg1, inarg2, wrkarg, outarg1, ... )
c   input argument list:
c     inarg1   - generic description, including content, units,
c     inarg2   - type.  explain function if control variable.
c
c   output argument list:      (including work arrays)
c     wrkarg   - generic description, etc., as above.
c     outarg1  - explain completely if error return
c     errflag  - even if many lines are needed
c
c   input files:   (delete if no input files in subprogram)
c     ddname1  - generic name & content
c
c   output files:  (delete if no output files in subprogram)
c     ddname2  - generic name & content as above
c     ft06f001 - include if any printout
c
c remarks: list caveats, other helpful hints or information
c
c attributes:
c   language: indicate extensions, compiler options
c   machine:  nas, cyber, whatever
c
c$$$
      subroutine gauss_lat_nmc(gaul,k)
c
      implicit real*8(a-h,o-z)
      dimension a(500)
      real*4 gaul(1)
c
      save
c
      esp=1.d-14
      c=(1.d0-(2.d0/3.14159265358979d0)**2)*0.25d0
      fk=k
      kk=k/2
      call bsslz1(a,kk)
      do 30 is=1,kk
      xz=cos(a(is)/sqrt((fk+0.5d0)**2+c))
      iter=0
   10 pkm2=1.d0
      pkm1=xz
      iter=iter+1
      if(iter.gt.10) go to 70
      do 20 n=2,k
      fn=n
      pk=((2.d0*fn-1.d0)*xz*pkm1-(fn-1.d0)*pkm2)/fn
      pkm2=pkm1
   20 pkm1=pk
      pkm1=pkm2
      pkmrk=(fk*(pkm1-xz*pk))/(1.d0-xz**2)
      sp=pk/pkmrk
      xz=xz-sp
      avsp=abs(sp)
      if(avsp.gt.esp) go to 10
      a(is)=xz
   30 continue
      if(k.eq.kk*2) go to 50
      a(kk+1)=0.d0
      pk=2.d0/fk**2
      do 40 n=2,k,2
      fn=n
   40 pk=pk*fn**2/(fn-1.d0)**2
   50 continue
      do 60 n=1,kk
      l=k+1-n
      a(l)=-a(n)
   60 continue
c
      radi=180./(4.*atan(1.))
      do 211 n=1,k
      gaul(n)=acos(a(n))*radi-90.0
  211 continue
c     print *,'gaussian lat (deg) for jmax=',k
c     print *,(gaul(n),n=1,k)
c
      return
   70 write(6,6000)
 6000 format(//5x,14herror in gauaw//)
      stop
      end

      subroutine gauss_lat_pcmdi(gaul,gaulb,k)
      parameter (nmax=500)
      real*8 pa(nmax),pw(nmax),pi,pb(nmax),pt(nmax),r2d
      real*4 gaul(k),gaulb(k)
C
C         get the gaussian latitudes and integration weights
C
CCC      call gauaw(pa,pw,k,nmax)
      call gauaw(pa,pw,k)
C
C         reverse direction so j increase northward for output
C
      do j=1,k
        jj=k-j+1
        pt(j)=pa(jj)
      end do
      do j=1,k
        pa(j)=pt(j)
      end do

      do j=1,k
        jj=k-j+1
        pt(j)=pw(jj)
      end do
      do j=1,k
        pw(j)=pt(j)
      end do

      pi = 4.*datan(1.d0)
      r2d = 180.0/pi
C
C         integrate to get the latitude boundaries of the gauss grid boxes
C
      pb(1)=-pi*0.5
      do j=1,k
        pb(j+1) = dasin( pw(j) + dsin( pb(j) ) )
      end do

      do j=1,k+1
        gaulb(j)=r2d*pb(j)
      end do

      do j=1,k
        gaul(j)=r2d*dasin(pa(j))
      end do

      return
      end

      subroutine gauaw(pa,pw,k)
c
c****     *gauaw* - compute abscissas and weights for *gaussian integration.
c
c         purpose.
c         --------
c
c         *gauaw* is called to compute the abscissas and weights requir
c         to perform *gaussian integration.
c
c**       interface.
c         ----------
c
c         *call* *gauaw(pa,pw,k)*
c
c         *pa*     - array, length at least *k,* to receive abscis
c         abscissas.
c         *pw*     - array, length at least *k,* to receive
c         weights.
c
c         method.
c         -------
c
c         the zeros of the *bessel functions are used as starting
c         approximations for the abscissas. newton iteration is used to
c         improve the values to within a tollerence of *eps.*
c
c         external.
c         ---------
c
c         *bsslzr* - routine to obtain zeros of *bessel functions.
c
c         reference.
c         ----------
c
      implicit real*8 (a-h, o-z)
      real*8 pa, pw
      dimension pa(k),pw(k)
      data eps/1.d-13/
c
c         ------------------------------------------------------------------
c
c*        1.     set constants and find zeros of bessel function.
c         --- --------- --- ---- ----- -- ------ ---------
c
      pi = 4.*datan(1.d0)
 100  continue
      c=(1.-(2./pi)**2)*0.25
      fk=k
      kk=k/2
      call bsslzr(pa,kk)
c
      do 290 is=1,kk
        xz=dcos(pa(is)/dsqrt((fk+0.5)**2+c))
c*        giving the first approximation for *xz.*
        iter=0
c
c         ------------------------------------------------------------------
c
c*        2.     compute abscissas and weights.
c         ------- --------- --- -------
c
 200    continue
c
c*        2.1     set values for next iteration.
 210    continue
        pkm2=1.
        pkm1=xz
        iter=iter+1
        if(iter.gt.10) go to 300
c
c*        2.2     computation of the *legendre polynomial.
 220    continue
c
        do 222 n=2,k
          fn=n
          pk=((2.*fn-1.)*xz*pkm1-(fn-1.)*pkm2)/fn
          pkm2=pkm1
          pkm1=pk
 222    continue
c
        pkm1=pkm2
        pkmrk=(fk*(pkm1-xz*pk))/(1.-xz**2)
        sp=pk/pkmrk
        xz=xz-sp
        avsp=dabs(sp)
        if(avsp.gt.eps) go to 210
c
c*        2.3     abscissas and weights.
 230    continue
        pa(is)=xz
        pw(is)=(2.*(1.-xz**2))/(fk*pkm1)**2
c
c*        2.4     odd *k* computation of weight at the equator.
 240    continue
        if (k.ne.kk*2) then
          pa(kk+1)=0.
          pk=2./fk**2
c
          do 242 n=2,k,2
            fn=n
            pk=pk*fn**2/(fn-1.)**2
 242      continue
c
          pw(kk+1)=pk
        else
c
c*        2.5     use symmetry to obtain remaining values.
c
 250      continue
c
          do 252 n=1,kk
            l=k+1-n
            pa(l)=-pa(n)
            pw(l)=pw(n)
 252      continue
c
        endif
 290  continue
c         write(6,*) iter
c
      return
c
c         ------------------------------------------------------------------
c
c*        3.     error processing.
c         ----- -----------
c
 300  continue
      write (nout,9901)
 9901 format(//,'  gauaw failed to converge after 10 iterations.')
      stop
c
c         ------------------------------------------------------------------
c
      end

      subroutine bsslzr(pbes,knum)
c
c****     *bsslzr* - routine to return zeros of the j0 *bessel function.
c
c         purpose.
c         --------
c
c         *bsslzr* returns *knum* zeros, or if *knum>50,* *knum*
c         approximate zeros of the *bessel function j0.
c
c**       interface.
c         ----------
c
c         *call* *nsslzr(pbes,knum)*
c
c         *pbes*   - array, dimensioned *knum,* to receive the
c         values.
c         *knum*   - number of zeros requested.
c
c         method.
c         -------
c
c         the first 50 values are obtained from a look up table. any
c         additional values requested are interpolated.
c
c         externals.
c         ----------
c
c         none.
c
c         reference.
c         ----------
c
*call     comcon
      real*8 pbes(knum), zbes(50), api
      data zbes        / 2.4048255577,   5.5200781103,
     x 8.6537279129,  11.7915344391,  14.9309177086,  18.0710639679,
     x 21.2116366299,  24.3524715308,  27.4934791320,  30.6346064684,
     x 33.7758202136,  36.9170983537,  40.0584257646,  43.1997917132,
     x 46.3411883717,  49.4826098974,  52.6240518411,  55.7655107550,
     x 58.9069839261,  62.0484691902,  65.1899648002,  68.3314693299,
     x 71.4729816036,  74.6145006437,  77.7560256304,  80.8975558711,
     x 84.0390907769,  87.1806298436,  90.3221726372,  93.4637187819,
     x 96.6052679510,  99.7468198587, 102.8883742542, 106.0299309165,
     x 109.1714896498, 112.3130502805, 115.4546126537, 118.5961766309,
     x 121.7377420880, 124.8793089132, 128.0208770059, 131.1624462752,
     x 134.3040166383, 137.4455880203, 140.5871603528, 143.7287335737,
     x 146.8703076258, 150.0118824570, 153.1534580192, 156.2950342685/
c
c         ------------------------------------------------------------------
c
c*        1.     extract values from look up table.
c         ------- ------ ---- ---- -- ------
c
c         set api
c
      api=4.*datan(1.d0)
 100  continue
      inum=min0(knum,50)
c
      do 110 j=1,inum
        pbes(j)=zbes(j)
 110  continue
c
c         ------------------------------------------------------------------
c
c*        2.     interpolate remaining values.
c         ----------- --------- -------
c
 200  continue
c
      zpi=api
      do 210 j=51,knum
        pbes(j)=pbes(j-1)+api
 210  continue
c
c         ------------------------------------------------------------------
c
      return
      end

      subroutine qprntu(a,qtitle,undef,
     $     ibeg,jbeg,m,n,iskip,iunit)
C
C----------------------------------------------------------------
C
C         version of qprntn which handles undefined values
C         by having them printed as ****
C
C         Mike Fiorino, NMC Development Division
C
C----------------------------------------------------------------
C
C
C         a= fwa of m x n array
C         qtitle - title
C         ibeg,jbeg=lower left corner coords to be printed
C         up to 43 x 83 points printed
c
      dimension a(m,n),ix(81)
      character qtitle*24
c
c  determine grid limits
c
      if(iskip.eq.0) iskip=1
      iend=min0(ibeg+79*iskip,m)
      jend=min0(jbeg+79*iskip,n)
c
   24 continue
c
c  index backwards checking for max
c
   11 xm=0.
      jendsc=min0(jend,n)
      do j=jbeg,jendsc,iskip
        jend_qp = j
        do i=ibeg,iend,iskip
          xmax=abs(a(i,j))
          if(a(i,j).eq.undef) xmax=0.0
          xm=amax1(xm,xmax)
        end do
      end do
c
c  determine scaling factor limits
c
      if(xm.lt.1.0e-32.or.xm.eq.0.0) xm=99.0
      xm=alog10(99.0/xm)
      kp=xm
      if(xm.lt.0.0)kp=kp-1
c
c  print scaling constants
c
   12 write(iunit,1) qtitle,kp,iskip,(i,i=ibeg,iend,2*iskip)

    1 format('0',a,'   k=',i3,' iskip=',i2,/,' ',41i6)
      fk=10.0**kp
c
c  quickprint field
c
      do 2 jli=jend_qp,jbeg,-iskip
        ii= 0
        if(kp.eq.0) then
          do i=ibeg,iend,iskip
            ii=ii+1

            if(a(i,jli).eq.undef) then
              ix(ii)=999999
            else
              ix(ii)=a(i,jli)+sign(.5,a(i,jli))
            endif

          end do
        else
          do i=ibeg,iend,iskip
            ii=ii+1
            if(a(i,jli).eq.undef) then
              ix(ii)=999999
            else
              ix(ii)=a(i,jli)*fk+sign(.5,a(i,jli))
            endif
          end do
        end if
        write(iunit,'(i4,81i3)') jli,(ix(i),i=1,ii),jli
2     continue
      return
      end

      subroutine qprntn(a,qtitle,ibeg,jbeg,m,n,iskip,iunit)
c
c**********	12 APR 91 this version outputs to iunit
c**********	using write on the Cray Y/MP
c
c***************************************************************
c***************************************************************
c*****                                                     *****
c*****       qprint output routine (corrected 4/26/86)     *****
c*****                                                     *****
c***************************************************************
c***************************************************************
c
c a= fwa of m x n array
c qtitle - title
c ibeg,jbeg=lower left corner coords to be printed
c up to 43 x 83 points printed
c
      dimension a(m,n),ix(81)
      character qtitle*24
c
c  determine grid limits
c
      if(iskip.eq.0) iskip=1
      iend=min0(ibeg+79*iskip,m)
      jend=min0(jbeg+79*iskip,n)
c
   24 continue
c
c  index backwards checking for max
c
   11 xm=0.
      jendsc=min0(jend,n)
      do j=jbeg,jendsc,iskip
      jend_qp = j
      do i=ibeg,iend,iskip
        xm=amax1(xm,abs(a(i,j)))
      end do
      end do
c
c  determine scaling factor limits
c
      if(xm.lt.1.0e-32.or.xm.eq.0.0) xm=99.0
      xm=alog10(99.0/xm)
      kp=xm
      if(xm.lt.0.0)kp=kp-1
c
c  print scaling constants
c
   12 write(iunit,1) qtitle,kp,iskip,(i,i=ibeg,iend,2*iskip)

    1 format('0',a,'   k=',i3,' iskip=',i2,/,' ',41i6)
      fk=10.0**kp
c
c  quickprint field
c
      do 2 jli=jend_qp,jbeg,-iskip
        ii= 0
        if(kp.eq.0) then
          do i=ibeg,iend,iskip
            ii=ii+1
            ix(ii)=a(i,jli)+sign(.5,a(i,jli))
          end do
        else
          do i=ibeg,iend,iskip
            ii=ii+1
            ix(ii)=a(i,jli)*fk+sign(.5,a(i,jli))
          end do
        end if
        write(iunit,'(i4,81i3)') jli,(ix(i),i=1,ii),jli
2     continue
      return
      end
      function qsatw(t,p)
c
c     t is temperature of air in deg celcius.
c     ta is temperature in deg kelvin
c     p is pressure in mb
c     qsatw is saturation specific humidity in g/g (over water)
c
      data ps/1013.246/,ts/373.16/
c
      ta = t
      if(t .lt. 100.0) then
         ta = t + 273.16
      end if
      e1=11.344*(1.0-ta/ts)
      e2=-3.49149*(ts/ta-1.0)
      f1=-7.90298*(ts/ta-1.0)
      f2=5.02808*alog10(ts/ta)
      f3=-1.3816*(10.0**e1-1.0)*1.e-7
      f4=8.1328*(10.0**e2-1.0)*1.e-3
      f5=alog10(ps)
      f=f1+f2+f3+f4+f5
      es=10.0**f
      qsatw=.62197*es/(p-0.378*es)
      return
      end
c
      function qsati(t,p)
c
c     t is temperature of air in deg celcius.
c     ta is temperature in deg kelvin
c     p is pressure in mb
c     qsati is saturation specific humidity with respect to ice in g/g
c
      data p0/6.1071/,t0/273.16/
c
      t = ta
      if(t .lt. 100.0) then
         ta = t + 273.16
      end if
      f1=-9.09718*(t0/ta-1.0)
      f2=-3.56654*alog10(t0/ta)
      f3=0.876793*(1.0-ta/t0)
      f4=alog10(p0)
      es=10.0**(f1+f2+f3+f4)
      qsati=.62197*es/(p-0.378*es)
      return
      end

      subroutine read_grid(iunit,a,ni,nj,istat)
      dimension a(ni,nj)
      istat=1
      read(iunit,err=800) a
      go to 999
 800  continue
      istat=0
 999  continue
      return
      end

      function satvp(temp)
      real a(7)
      data a/6.10779961,0.4436518521,0.01428945805,2.650648471e-4,
     &       3.031240396e-6,2.034080948e-8,6.136820929e-11/
      data e0, aa, bb/6.1078, 17.2693882, 237.3/
      data t0/273.16/
c
      t = temp
c convert temperature to degrees c, if necessary.
      if(temp .gt. 100.0) t = temp - t0
c polynomial is only good from -50 to + 50 c.
c if outside this range teton's formmula is used.
      if(t .gt. 50.0 .or. t .lt. -50.0) then
        satvp = e0*exp(aa*t/(t+bb))
      else
      satvp = a(1)+t*(a(2)+t*(a(3)+t*(a(4)+t*(a(5)+t*(a(6)+a(7)*t)))))
      end if
      return
      end

      subroutine sfc_area(fld,rlon,rlat,undef,ni,nj,
     $     area,iunit_diag)
      dimension fld(ni,nj),rlat(nj+1),rlon(ni+1),
     $     area(ni,nj)
      deg2rad=3.14115926/180.0
      do j=1,nj
        do i=1,ni
          dlon=(rlon(i+1)-rlon(i))*deg2rad
          dlat=sin(rlat(j+1)*deg2rad)-sin(rlat(j)*deg2rad)
          area(i,j)=dlon*dlat
          if(fld(i,j).eq.undef) area(i,j)=0.0
        end do
        rlatout=(rlat(j+1)+rlat(j))*0.5
cccc        write(iunit_diag,*) 'j = ',j,' sfc area = ',rlatout,area(1,j)
      end do

      return
      end

      subroutine trim_grid(a,ni,niinew,nj,dum)
      dimension a(ni,nj),dum(niinew,nj)
      do j=1,nj
        do i=1,niinew
          dum(i,j)=a(i,j)
        end do
      end do
      return
      end

      subroutine vload2(a,b,ni,nj)
      dimension a(ni,nj),b(ni,nj)
      do j=1,nj
        do i=1,ni
          a(i,j)=b(i,j)
        end do
      end do
      return
      end
      function ichar_len(c,imax)
      character*1 c(imax)
      iend=-1
      ii=1
      do while (iend.eq.-1.and.ii.le.imax)
        if(c(ii).eq.' ') iend=ii
        ii=ii+1
      end do
      if(ii.gt.imax) then
        ichar_len=imax
      else
        ichar_len=iend-1
      end if
      return
      end

      subroutine ul_case(cc,ilen)
      character*1 cc(ilen)
      do i=1,ilen
        if(ichar(cc(i)).ge.65.and.ichar(cc(i)).le.90)
     $       cc(i)=char(ichar(cc(i))+32)
      end do
      return
      end

      subroutine write_grid(iunit,a,ni,nj)
      dimension a(ni,nj)
      write(iunit) a
      return
      end

ext: f


text:
!===============================================================================
!  S H E B A - Shear-wave Birefringence Analysis
!===============================================================================
!  This software is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
!===============================================================================
!
!  James Wookey, School of Earth Sciences, University of Bristol
!

!=======================================================================
      subroutine read_config()
!=======================================================================
!
!     Read SHEBA configuration files from input file
!
      use sheba_config
      implicit none

      character :: dum ! just a dummy character
      character (len=80) :: fn

!  ** open file
      open(10,file='sheba.in',status='old',err=900)
!  ** skip header line
      read(10,*) dum

!  ** read input parameters
      read(10,*,err=901) fn
      config % fname_base = trim(fn)
      read(10,*,err=901) config % comp(1)
      read(10,*,err=901) config % comp(2)
      read(10,*,err=901) config % comp(3)
      read(10,*,err=901) config % imode

      if (config % imode == 0) then
         read(10,*,err=901) config % source_pol
      else
         config % source_pol = 0.0
      endif
      read(10,*,err=901) config % nwindbeg, config % nwindend
      read(10,*,err=901) config % max_tlag

!  ** read REC correction options
      read(10,*,err=901) config % iuma_corr
      if (config % iuma_corr == 1) then
         read(10,*,err=901) config % uma_phi,config % uma_dt
      else
         config % uma_phi = 0.0 ; config % uma_dt = 0.0
      endif

!  ** read SRC correction options
      read(10,*,err=901) config % i_src_corr
      if (config % i_src_corr == 1) then
         read(10,*,err=901) config % src_fast, config % src_tlag
      else
         config % src_fast = 0.0 ; config % src_tlag = 0.0
      endif

      return
!  ** error handling
900   print*,'<!> ERROR IN SHEBA: sheba.in not found'
      stop
901   print*,'<!> ERROR IN SHEBA: error reading sheba.in'
      stop
      end subroutine read_config
!=======================================================================

!=======================================================================
      subroutine get_traces(h1,h2,v)
!=======================================================================
!
!     Read the input traces and determine the h1,h2,v ordering
!
!-----------------------------------------------------------------------
      use f90sac  ! use the f90sac module
      use sheba_config ! use the sheba_config module
      use event_info ! use the event info module
!-----------------------------------------------------------------------
      implicit none
      type (SACTrace) :: h1,h2,v ! the re-ordered traces for analysis
      type (SACTrace) :: traces(3) ! the input traces
      character (len=80) :: fn ! filename
!      integer :: iorder(3) ! index to re-order the traces
!                 now in config
      integer :: i,ii

!  ** read the three files
      do i = 1,3
         fn = trim(config % fname_base) // '.' // config % comp(i)
         call f90sac_readtrace(fn,traces(i))
      enddo ! do i = 1,3

!  ** check that cmpaz and cmpinc are set
      if ( traces(1)%cmpaz == SAC_rnull .or.  &
           traces(2)%cmpaz == SAC_rnull .or.  &
           traces(3)%cmpaz == SAC_rnull) then
         print*,'<!> ERROR IN SHEBA: cmpaz header not set in one or more files'
         STOP
      endif
      if ( traces(1)%cmpinc == SAC_rnull .or.  &
           traces(2)%cmpinc == SAC_rnull .or.  &
           traces(3)%cmpinc == SAC_rnull) then
         print*,'<!> ERROR IN SHEBA: cmpinc header not set in one or more files'
         STOP
      endif

!  ** check that traces are the same length, delta and origin time
      if ( traces(1)%npts /= traces(2)%npts .or.  &
           traces(1)%npts /= traces(3)%npts ) then
         print*,'<!> ERROR IN SHEBA: Traces are different lengths'
         STOP
      endif
      if ( traces(1)%delta /= traces(2)%delta .or.  &
           traces(1)%delta /= traces(3)%delta ) then
         print*,'<!> ERROR IN SHEBA: Traces have different sampling rates'
         STOP
      endif

      if ( f90sac_compare_origin_time(traces(1),traces(2)) /= 0 .or.  &
           f90sac_compare_origin_time(traces(1),traces(3)) /= 0 ) then
         print*,'<!> ERROR IN SHEBA: Traces have different origin times'
         print*,traces(1)%nzyear,traces(1)%nzjday,traces(1)%nzhour,&
                traces(1)%nzmin,traces(1)%nzsec,traces(1)%nzmsec
         print*,traces(2)%nzyear,traces(2)%nzjday,traces(2)%nzhour,&
                traces(2)%nzmin,traces(2)%nzsec,traces(2)%nzmsec
         print*,traces(3)%nzyear,traces(3)%nzjday,traces(3)%nzhour,&
                traces(3)%nzmin,traces(3)%nzsec,traces(3)%nzmsec
         STOP
      endif

!  ** determine ordering and check inclinations
      ii = 1
      do i = 1,3
!     ** check component inclinations are 0 or 90

         if ( abs(traces(i) % cmpinc) > angtol .and.  &
             abs(traces(i) % cmpinc - 90.0) > angtol ) then
            print*,'<!> ERROR IN SHEBA: input trace is not ', &
               'horizontal or vertical'
            stop
         endif

         if (abs(traces(i) % cmpinc) < angtol) then
            config % iorder(3) = i
         else
            config % iorder(ii) = i
            ii = ii + 1
         endif
      enddo ! do i = 1,3

!  ** make h1 the smallest azimuth
      if (traces(config % iorder(1)) % cmpaz > traces(config % iorder(2)) % cmpaz) then
         i=config % iorder(1)
         config % iorder(1) = config % iorder(2)
         config % iorder(2) = i
      endif

!  ** assign the traces
      h1 = traces(config % iorder(1))
      h2 = traces(config % iorder(2))
      v = traces(config % iorder(3))

!  ** get some useful event information from the vertical header
      event % dist = v % dist
      event % az = v % az
      event % baz = v % baz

!  ** check that the horizontal traces are orthogonal
      if (f90sac_orth2d(h1,h2) == 0) then
         print*,'<!> ERROR IN SHEBA: input horizontals are ', &
               'not orthogonal'
         STOP
      endif
!  ** done!
      return
      end subroutine get_traces
!=======================================================================

max_stars_repo_name: jwookey/sheba

max_issues_repo_name: jwookey/sheba

max_forks_repo_name: jwookey/sheba


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

   SUBROUTINE libxc_postprocess_gga_mt(xcpot,atoms,sym,sphhar,noco,n,v_xc,grad, atom_num)
      USE m_mt_tofrom_grid
      USE m_types
      use m_judft_string

      IMPLICIT NONE
      CLASS(t_xcpot),INTENT(IN)   :: xcpot
      TYPE(t_atoms),INTENT(IN)    :: atoms
      TYPE(t_sym),INTENT(IN)      :: sym
      TYPE(t_sphhar),INTENT(IN)   :: sphhar
      TYPE(t_noco),INTENT(IN)     :: noco
      INTEGER,INTENT(IN)          :: n
      REAL,INTENT(INOUT)          :: v_xc(:,:)
      TYPE(t_gradients),INTENT(IN):: grad
      INTEGER, OPTIONAL           :: atom_num

      INTEGER :: nsp,n_sigma,i
      REAL,ALLOCATABLE:: vsigma(:,:),vsigma_mt(:,:,:)
      TYPE(t_gradients)::grad_vsigma
      character(len=:), allocatable :: fname

      n_sigma=MERGE(1,3,SIZE(v_xc,2)==1) !Number of contracted gradients in libxc 1 for non-spin-polarized, 3 otherwise
      nsp=SIZE(v_xc,1) !no of points
      ALLOCATE(vsigma(nsp,n_sigma),vsigma_mt(atoms%jri(n),0:sphhar%nlhd,n_sigma))
      vsigma_mt=0.0
      vsigma=TRANSPOSE(grad%vsigma) !create a (nsp,n_sigma) matrix
      CALL mt_from_grid(atoms,sym,sphhar,n,n_sigma,vsigma,vsigma_mt)
      DO i=1,atoms%jri(n)
         vsigma_mt(i,:,:)=vsigma_mt(i,:,:)*atoms%rmsh(i,n)**2
      ENDDO
      ALLOCATE(grad_vsigma%gr(3,nsp,n_sigma))
      CALL mt_to_grid(xcpot%needs_grad(),n_sigma,atoms,sym,sphhar,.TRUE.,vsigma_mt,n,noco,grad=grad_vsigma)

      CALL libxc_postprocess_gga(transpose(grad%vsigma),grad,grad_vsigma,v_xc)
   END SUBROUTINE libxc_postprocess_gga_mt

   SUBROUTINE libxc_postprocess_gga_pw(xcpot,stars,cell,v_xc,grad)
      USE m_pw_tofrom_grid
      USE m_types

      IMPLICIT NONE
      CLASS(t_xcpot),INTENT(IN)   :: xcpot
      TYPE(t_stars),INTENT(IN)    :: stars
      TYPE(t_cell),INTENT(IN)     :: cell
      REAL,INTENT(INOUT)          :: v_xc(:,:)
      TYPE(t_gradients),INTENT(IN):: grad

      COMPLEX,ALLOCATABLE:: vsigma_g(:,:)
      REAL,ALLOCATABLE:: vsigma(:,:)
      TYPE(t_gradients)::grad_vsigma
      INTEGER :: nsp,n_sigma

      nsp=SIZE(v_xc,1) !no of points
      n_sigma=MERGE(1,3,SIZE(v_xc,2)==1) !See in _mt routine
      ALLOCATE(vsigma_g(stars%ng3,n_sigma),vsigma(nsp,n_sigma)); vsigma_g=0.0
      vsigma=TRANSPOSE(grad%vsigma) !create a (nsp,n_sigma) matrix
      CALL pw_from_grid(xcpot%needs_grad(),stars,.FALSE.,vsigma,vsigma_g)
      !vsigma_g(:,1)=vsigma_g(:,1)*stars%nstr(:)
      ALLOCATE(grad_vsigma%gr(3,nsp,n_sigma))
      CALL pw_to_grid(xcpot%needs_grad(),n_sigma,.false.,stars,cell,vsigma_g,grad_vsigma,xcpot)

      CALL libxc_postprocess_gga(transpose(grad%vsigma),grad,grad_vsigma,v_xc)
   END SUBROUTINE libxc_postprocess_gga_pw

   SUBROUTINE libxc_postprocess_gga(vsigma,grad,grad_vsigma,v_xc)
      USE m_types
      IMPLICIT NONE
      REAL,INTENT(IN)             :: vsigma(:,:)
      TYPE(t_gradients),INTENT(IN):: grad,grad_vsigma
      REAL,INTENT(INOUT)          :: v_xc(:,:)
      INTEGER:: i
      IF (SIZE(v_xc,2)==1) THEN !Single spin
         DO i=1,SIZE(v_xc,1) !loop over points
            v_xc(i,1)=v_xc(i,1)-2*dot_PRODUCT(grad_vsigma%gr(:,i,1),grad%gr(:,i,1))-2*vsigma(i,1)*grad%laplace(i,1)
         ENDDO
      ELSE  !two spins
         DO i=1,SIZE(v_xc,1) !loop over points
            v_xc(i,1)=v_xc(i,1)-2*dot_PRODUCT(grad_vsigma%gr(:,i,1),grad%gr(:,i,1))-2*vsigma(i,1)*grad%laplace(i,1)-&
                       dot_PRODUCT(grad_vsigma%gr(:,i,2),grad%gr(:,i,2))-vsigma(i,2)*grad%laplace(i,2)
            v_xc(i,2)=v_xc(i,2)-2*dot_PRODUCT(grad_vsigma%gr(:,i,3),grad%gr(:,i,2))-2*vsigma(i,3)*grad%laplace(i,2)-&
                       dot_PRODUCT(grad_vsigma%gr(:,i,2),grad%gr(:,i,1))-vsigma(i,2)*grad%laplace(i,1)
         ENDDO
      END IF

   END SUBROUTINE libxc_postprocess_gga

END MODULE m_libxc_postprocess_gga

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text:
      subroutine genadj (colm,         rowp, icnt )
c
      use pointer_data
c
      include "common.h"
c
      integer rowp(nshg*nnz),         colm(nshg+1)
      integer adjcnt(nshg),    row_fill_list(nshg,15*nnz), mloc(1)
c                                          change ^ if overflow
c                                   also change overflow check in asadj TWICE
      integer tmprdim(1), nnonzero
      real*8, allocatable, dimension(:) :: tmpr

      adjcnt=0

      do iblk = 1, nelblk
c
c.... set up the parameters
c
         iel    = lcblk(1,iblk)
         lelCat = lcblk(2,iblk)
         lcsyst = lcblk(3,iblk)
         iorder = lcblk(4,iblk)
         nenl   = lcblk(5,iblk) ! no. of vertices per element
         nshl   = lcblk(10,iblk)
         npro   = lcblk(1,iblk+1) - iel

c
c.... compute sparse matrix data structures
c
         call Asadj (row_fill_list,
     &               mien(iblk)%p,  adjcnt )

      enddo

      call sumgatInt ( adjcnt, nshg, nnonzero)
      if ( myrank .eq. master) then
         write (*,*) 'Number of global nonzeros ',nnonzero
      endif

c
c     build the colm array
c
      colm(1)=1
      do i=1,nshg
         colm(i+1)=colm(i)+adjcnt(i)
      enddo
c
c     sort the rowp into increasing order
c
      ibig=10*nshg
      icnt=0
      tmprdim=maxval(adjcnt)
      allocate (tmpr(tmprdim(1)))
      do i=1,nshg
         ncol=adjcnt(i)
         tmpr(1:ncol)=row_fill_list(i,1:ncol)
         do j=1,ncol
            icnt=icnt+1
            imin=minval(tmpr(1:ncol))
            mloc=minloc(tmpr(1:ncol))
            rowp(icnt)=imin
            tmpr(mloc(1))=ibig
         enddo
      enddo
      deallocate(tmpr)
c      maxfill=tmprdim(1)
c      write(*,*) 'maxfill=',maxfill
      nnza=icnt/nshg +1
      if(icnt.gt.nnz*nshg) then
         write(*,*) 'increase nnz in genmat to',nnza
         stop
c      else
c         write(*,*) 'nnz ok  nnz=',nnz,' actually needed',nnza
c         write(*,*) myrank,' is my rank and my nnz_tot is: ',icnt
      endif
      return
      end

max_stars_repo_name: polmes/phasta

max_issues_repo_name: polmes/phasta

max_forks_repo_name: polmes/phasta


text:
! Copyright (C) 2015-2017, Simone Marsili
! All rights reserved.
! License: BSD 3 clause

module fasta
  use constants
  implicit none
  private
  public :: fasta_read
  public :: set_fasta_alphabet
  public :: fasta_alphabet
  integer, parameter         :: kprot=21, knuc=6
  ! TODO: "gap" symbol should map to index 1
  character(len=kprot)       :: protein_alphabet='-ACDEFGHIKLMNPQRSTVWY'
  character(len=knuc)        :: nuc_acid_alphabet='-ACGTU'
  character(len=string_size) :: fasta_alphabet

contains

  subroutine set_fasta_alphabet(data_type)
    character(len=*), intent(in) :: data_type

    select case(data_type)
    case('protein')
       fasta_alphabet = protein_alphabet
    case('nuc_acid')
       fasta_alphabet = nuc_acid_alphabet
    case default
       write(0,*) 'warning: not a valid alphabet for bioseqs'
    end select

  end subroutine set_fasta_alphabet

  function is_nuc_acid(string) result(res)
    ! check if a string contains only natural nucleotides symbols
    character(len=*), intent(in) :: string
    logical :: res

    if (verify(trim(string), trim(nuc_acid_alphabet)) == 0) then
       res = .true.
    else
       res = .false.
    end if

  end function is_nuc_acid

  function is_protein(string) result(res)
    ! check if a string contain only natural amino acids symbols
    character(len=*), intent(in) :: string
    logical :: res

    if (verify(trim(string), trim(protein_alphabet)) == 0) then
       res = .true.
    else
       res = .false.
    end if

  end function is_protein

  subroutine fasta_string2seq(alphabet, string,seq)
    ! read a sequence, return an array of integer
    character(len=*), intent(in)  :: alphabet
    character(len=*), intent(in)  :: string
    integer,          intent(out) :: seq(:)
    integer :: i

    seq = [(index(trim(alphabet), string(i:i)), i = 1, len_trim(string))]

  end subroutine fasta_string2seq

  subroutine fasta_read(unt, data, data_type, error_code)
    ! read n sequences from unit udata
    integer,              intent(in)  :: unt
    integer, allocatable, intent(out) :: data(:, :)
    character(len=*),     intent(out) :: data_type
    integer,              intent(out) :: error_code
    character(len=long_string_size) :: header
    character(len=long_string_size) :: string
    character(len=long_string_size) :: line
    integer, allocatable            :: seq(:)
    integer                         :: k, ns, nv, nnuc, nprot, ntot, nuc, prot
    integer                         :: err

    ! count data and Ls
    ns = 0
    ntot = 0
    nnuc = 0
    nprot = 0
    nuc = 0
    prot = 0
    do
       read(unt,'(a)',iostat=err) line
       ! exit at end of file or empty line
       if ( err < 0 .or. len_trim(line) == 0 ) exit
       if ( line(1:1) == ">" ) then
          ntot = ntot + 1
          if (ntot > 1) then
             nnuc = nnuc + nuc
             nprot = nprot + prot
          end if
       else
          if (is_nuc_acid(line)) nuc = 1
          if (is_protein(line)) prot = 1
       end if
    end do
    rewind(unt)
    nnuc = nnuc + nuc
    nprot = nprot + prot

    ! set data type
    if (nnuc > nprot) then
       data_type = 'nuc_acid'
       ns = nnuc
    else
       data_type = 'protein'
       ns = nprot
    end if

    call set_fasta_alphabet(data_type)

    string = ""
    k = 1
    read(unt,'(a)',iostat=err) header
    do
       read(unt,'(a)',iostat=err) line
       if (err < 0) then
          ! eof
          if (k == 1) then
             write(0,*) 'ERROR ! empty alignment'
             error_code = 1
          end if
          call fasta_string2seq(fasta_alphabet,string,seq)
          if (all(seq > 0)) data(:,k) = seq
          exit
       end if
       if (line(1:1) == ">") then
          header = line
          if (header(1:1) /= '>') then
             write(0,*) 'ERROR ! not a valid FASTA file'
             error_code = 2
          end if
          if (k == 1) then
             nv = len_trim(string)
             allocate(seq(nv), stat=err)
             allocate(data(nv, ns), stat=err)
          end if
          call fasta_string2seq(fasta_alphabet, string,seq)
          string = ""
          if (all(seq > 0)) then
             data(:, k) = seq
             k = k + 1
          end if
          cycle
       else
          string = trim(string)//trim(line)
       end if
    end do

  end subroutine fasta_read

end module fasta

max_stars_repo_name: simomarsili/elss

max_issues_repo_name: simomarsili/elss

max_forks_repo_name: simomarsili/elss