text:
C     LAST UPDATE 16/03/89
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      SUBROUTINE GETINS (ITERM,IPRINT,IOPT)
      IMPLICIT NONE
C
C Purpose: Read instruction from terminal, check to make sure its in
C          upper case and then calculate its option number from its
C          position in the command string.
C
      INTEGER ITERM,IPRINT,IOPT
C
C ITERM  : Terminal input stream
C IPRINT : Terminal output stream
C IOPT   : Instruction option number 0 - <ctrl-z>
C                                   >1 - valid instruction
C
C Calls   2: ERRMSG , UPCASE
C Called by: BSL
C
C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C LOCAL VARIABLES:
C
      INTEGER       INDX
      CHARACTER*256 CMDSTR
      CHARACTER*4   TBUFF
      LOGICAL       IGETS
C
C INDX   : Postion of valid instruction in command string
C CMDSTR : Command string
C TBUFF  : Terminal buffer
C
      DATA CMDSTR(1:36)    /'.ADC.ADD.ADN.ASF.AVE.CON.CUT.DIC.DIN'/
      DATA CMDSTR(37:72)   /'.DIS.DIV.DSF.DUP.EXP.FFT.HOR.IFT.INT'/
      DATA CMDSTR(73:108)  /'.IPO.ITP.LOG.MAX.MIR.MSF.MSK.MUC.MUL'/
      DATA CMDSTR(109:144) /'.MUN.PAK.POL.POW.PRT.CIR.REM.REP.RMP'/
      DATA CMDSTR(145:180) /'.ROT.SHF.SUM.SUN.VER.WIN.ZER.BAK.CHG'/
      DATA CMDSTR(181:216) /'.DGL.ARG.LZR.SEC.GAU.CNV.MRG.INS.ROI'/
      DATA CMDSTR(217:252) /'.CIN.RIN.SUR.RAD.AZI.***.***.***.***'/
      DATA CMDSTR(253:256) /'.***'/
C
C-----------------------------------------------------------------------
10    WRITE (IPRINT,'(A,$)') ' Enter instruction: '
      IF (IGETS (ITERM,TBUFF)) THEN
         CALL UPCASE (TBUFF,TBUFF)
         INDX=INDEX (CMDSTR,TBUFF)
         IF (INDX.NE.0) THEN
            IOPT=(INDX+3)/4
            RETURN
         ENDIF
         CALL ERRMSG ('Error: Illegal instruction')
         GOTO 10
      ENDIF
      IOPT=0
C
      END


text:
C    @(#)bldtbx.f	20.7 1/4/99
        subroutine bldtbx (nb, jtb, jxb, kerr)
c
c       This routine builds a complete entry in TBX corresponding to
c       data for bus "NB".
c       input data:
c
c               nb:   Bus number of completed and linked "busdta" and
c                     "bctbl",and if applicable,"xdata" arrays
c               jtb:  Tbx index in which entry will be made
c               jxb:  Xdata index of pertinent xdata (if applicable)
c
c       output data:
c
c               tbx:  Complete data assigned
c               kerr: Error status switch
c                       0 - normal
c                       1 - minor errors
c                       2 - fatal error
c
      include 'ipfinc/parametr.inc'

      include 'ipfinc/alpha.inc'
      include 'ipfinc/blank.inc'
      include 'ipfinc/bus.inc'
      include 'ipfinc/cbus.inc'
      include 'ipfinc/com008.inc'
      include 'ipfinc/lfiles.inc'
      include 'ipfinc/prt.inc'
      include 'ipfinc/tbx.inc'
      include 'ipfinc/xdata.inc'

      common /dbkk/ dgkk, dbkk, bkkadj
      double precision dgkk,dbkk, bkkadj

c     Local variables changed to double precision

      double precision qgen, qmax, qmin, qadj, totcap, totrek,
     &       usecap, userek
      character xbuf*120
      integer oldltp, olditp, bustype, status, ptr
c
c       build TBX entries according to type
c
        bustype = kbsdta(1,nb)
        ltype = tbxtyp(bustype)
        itype = 1
        if (ltype .eq. 0) then
           write (errbuf(1),130)
  130      format(' Illegal bus type call to "BLDTBX".')
           errbuf(2) = ' '
           call bcdbus (nb,xbuf)
           write (errbuf(3),150) xbuf(1:80)
  150      format(2x,a80)
           call prterx ('W',3)
           go to 230
        endif
c
c       search for entry
c
        ntotb = ntotb + 1
        if (ntotb .gt. MAXTBX) then
           write (errbuf(1),60) MAXTBX,bus(nb),base(nb)
   60      format (' More than ',i4,' special buses (types ',
     &             '"BV", "BQ", "BG", "BO", "BX", "BF".  ',
     &             'Overflow encountered at bus ',a8,f6.1)
           call prterx ('W',1)
           ntotb=1
        endif
        jtb=ntotb
c
c       obtain bus quantities
c
        qmax = busdta(9,nb)/bmva
        qmin = busdta(10,nb)/bmva
        ptr = kbsdta(15,nb)
        do while (ptr .gt. 0)
          qmax = qmax + bctbl(11,ptr)/bmva
          qmin = qmin + bctbl(12,ptr)/bmva
          ptr = bctbl_nxt(ptr)
        enddo

        qgen = qnetu(nb) + qloadu(nb)

        go to(170,180,190,180,190,170) ltype
c
c       "BV" and "BF"
c
  170   tbx(8,jtb)=0

        tbx(3,jtb)=vlimx(nb)
        tbx(4,jtb)=vlimn(nb)
        tbx(5,jtb)=(qgen-qloadu(nb))
        tbx(6,jtb)=0.0
        go to 206
c
c       "BQ", "BO", and "BL"
c
  180   tbx(8,jtb)=0
        qadj=bkkadj*vstart(nb)**2
        tbx(3,jtb)=qmax-qloadu(nb)-dmin1(0.0d0, qadj)
        tbx(4,jtb)=qmin-qloadu(nb)-dmax1(0.0d0, qadj)
        tbx(5,jtb)=vstart(nb)
        tbx(6,jtb)=bkkadj
        go to 202
c
c       "BG" and "BX" bus
c
  190   tbx(8,jtb)=kbsdta(13,nb)
        tbx(3,jtb)=qmax-qloadu(nb)
        tbx(4,jtb)=qmin-qloadu(nb)
        if (ltype .eq. 3) then
          tbx(5,jtb)=busdta(14,nb)
c
c         Reset % to zero if qmin=qmax
c
          if (tbx(3,jtb) - tbx(4,jtb) .lt. 0.01) then
            write (errbuf(1),10203) qmax,qmin
10203       format('Insufficent var range: QMAX (', f6.1,
     &             '), QMIN (', f6.1,') is <= 1 MVAR.')
            write (errbuf(2),10204)
10204       format('% is set to 0.0 on the following bus')
            call bcdbus (nb,xbuf)
            write (errbuf(3),150) xbuf(1:80)
            call prterx ('I',3)
            tbx(5,jtb) = 0.0
          endif
          tbx(6,jtb)=bkkadj
        else
          tbx(5,jtb)=jxb
          tbx(6,jtb)=bkkadj
c
c         Type X-bus may be initialized in a state other than ITYPE = 1.
c
          totrek = xdata(3,jxb)
          totcap = xdata(4,jxb)
          userek = xdata(5,jxb)
          usecap = xdata(6,jxb)
          if (totrek .ne. 0.0 .and. totcap .ne. 0.0) then
            if (totrek .eq. userek .and. totcap .eq. userek) then
              itype = 1
            else if (userek .eq. 0.0 .and. usecap .eq. totcap) then
              itype = 3
            else if (usecap .eq. 0.0 .and. userek .eq. totrek) then
              itype = 2
            else
              itype = 4
            endif
          else if (totrek .ne. 0.0) then
            if (userek .eq. totrek) then
              itype = 2
            else if (userek .eq. 0.0) then
              itype = 3
            else
              itype = 4
            endif
          else
            if (usecap .eq. totcap) then
              itype = 3
            else if (usecap .eq. 0.0) then
              itype = 2
            else
              itype = 4
            endif
          endif
        endif
c
c       check q limits
c
  202   if (qmax .lt. qmin) then
           if (abs (qmax-qmin) .lt. 0.1) then
              busdta(10,nb) = qmax
           else
              write (errbuf(1),204) qmax,qmin
  204         format('0 "qmax" (',f6.1,') is less than "qmin" (',
     1                f6.1,') ', 'on the following bus:')
              errbuf(2) = ' '
              call bcdbus (nb,xbuf)
              write (errbuf(3),150) xbuf(1:80)
              call prterx ('W',3)
           endif
        endif
c
c       Reset "TBX" for hot start
c
  206   call getoldtb (bus(nb), base(nb), oldjtb, oldltp, olditp,
     &                 status)
        if (status .eq. 1) then
           if (ltype .eq. oldltp .and. itype .ne. olditp) then
              itype = olditp
              call updatetb(nb, oldjtb, jtb)
           endif
        endif
c
c       complete "tbx"
c
        tbx(1,jtb) = ltype
        tbx(2,jtb) = nb
        tbx(7,jtb) = itype

  230   return
        end

ext: f

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
      subroutine diag_memory_all(n_cf,ncoef,ico,nz,hmx,
     :           coeff,value,inptr,jptr,jjh,nijcurr,nze)

      implicit double precision(a-h,o-z)
      parameter(MEIG=20,MTERM=20)

      dimension ico(nze), hmx(nze), coeff(*), value(*), inptr(*)
      dimension jptr(*)
      hmx(1:nze) = 0.0

      max_col = 1;
      do ii = 1, n_cf;
*       .. test for next non-zero matrix element
          n_count_tmp = ncoef+ii
          if (ii.gt.ico(nijcurr+nz)) nijcurr = nijcurr + 1
          hmx(nijcurr) = hmx(nijcurr) +
     :            coeff(n_count_tmp)*value(inptr(n_count_tmp))
          if (nijcurr.gt.jptr(jjh)) then
            jjh = jjh + 1;
            max_col = max_col+1
          end if
*          print '(3I5,3F20.15,I4)',
*     :     ii,nijcurr,inptr(ncoef+ii),coeff(ncoef+ii),
*     :                value(inptr(ncoef+ii)),
*     :     hmx(nijcurr)
        end do

      ncoef = ncoef + n_cf;

      end

****
      subroutine diag_disk_clst(ncoef,ico,hmx,coeff,value,
     :           inptr,nijcurr,ncodim,nze,nz)

      implicit double precision(a-h,o-z)
      parameter(MEIG=20,MTERM=20)

      dimension ico(nz), hmx(nz), coeff(*), value(*), inptr(*)

      hmx(1:nze) = 0
      ncoef = 0;
        do
          read(39) num,(coeff(j),j=1,num),
     :                 (inptr(j),j=1,num)
          do ii = 1, num
*       .. test for next non-zero matrix element
            if (ncoef+ii.gt.ico(nijcurr+nz)) nijcurr = nijcurr + 1
            hmx(nijcurr) = hmx(nijcurr) + coeff(ii)*value(inptr(ii))
*           print '(3I5,3F20.15,I4)',
*     :     ii,nijcurr,inptr(ii),coeff(ii),value(inptr(ii)),
*     :     hmx(nijcurr)
          end do
          ncoef = ncoef + num
          if (num .lt. ncodim) exit
        end do
*     .. compute diagonals

        end

*****
      subroutine diag_disk_ico(ncoef,ico,hmx,coeff,value,inptr,
     :           jptr,jjh,nijcurr,ncodim,ncfg,nze,nze_col)

      implicit double precision(a-h,o-z)
      parameter(MEIG=20,MTERM=20)
      logical write_hmx
      dimension ico(nze_col),hmx(nze), coeff(*), value(*), inptr(*)
      dimension jptr(ncfg);

      ncoef = 0;
      nijcurr = 1;
      nij_count = 1;
      hmx(1:nze) = 0.0
      read(12) nico, (ico(jj),jj=1,nico)
      itmp = jptr(1)
      write_hmx = .false.
      rewind(13)

      do
         read(39) num,(coeff(j),j=1,num), (inptr(j),j=1,num)

         do ii = 1, num

           if (ncoef+ii > ico(nij_count)) then
              nijcurr = nijcurr + 1
              nij_count = nij_count + 1
           end if

           hmx(nijcurr) = hmx(nijcurr) +
     :                        coeff(ii)*value(inptr(ii))

           if (ncoef+ii > ico(nico)) then
              if (jjh < ncfg) then
                 read(12) nico ,(ico(jj),jj=1,nico)
                 nij_count = 1
              end if
              if (nijcurr-1 == jptr(jjh)) then
                 jjh = jjh + 1
              end if
           end if

*           print '(3I5,3F20.15,I4)',
*     :     ii,nijcurr,inptr(ii),coeff(ii),value(inptr(ii)),
*     :     hmx(nijcurr)
         end do
         ncoef = ncoef + num
         if (num .lt. ncodim) exit
      end do
      End subroutine diag_disk_ico

*************
      subroutine diag_disk_hmx(ncoef,ico,hmx,coeff,value,inptr,
     :           jptr,jjh,nijcurr,ncodim,ncfg,nze_max,nze,hii,
     :           shift_hmx)

      implicit double precision(a-h,o-z)
      parameter(MEIG=20,MTERM=20)
      logical next_nze,read_ico
      dimension ico(nze_max), hmx(nze_max), coeff(*), value(*), inptr(*)
      dimension jptr(ncfg), hii(ncfg);

      ncoef = 0;
      nijcurr = 1;
      nij_count = 1;
      hmx(1:nze_max) = 0.0
      ico(1:nze_max) = 0
      read(12) nico, (ico(jj),jj=1,nico)
      itmp = jptr(1)
      read_ico = .false.
      next_nze = .false.
      shift_hmx = 0
      rewind(13)

      do
         read(39) num,(coeff(j),j=1,num), (inptr(j),j=1,num)

         do ii = 1, num
           hmx(nij_count) = hmx(nij_count) +
     :                        coeff(ii)*value(inptr(ii))
           if (ncoef+ii == ico(nij_count))  next_nze = .true.

*           print '(3I5,3F20.15,I4)',
*     :     ii,nij_count,inptr(ii),coeff(ii),value(inptr(ii)),
*     :     hmx(nij_count)

           if ((next_nze).and.(nij_count == 1)) then
              if (jjh == 1) shift_hmx = hmx(1)
              hii(jjh)= hmx(nij_count) - shift_hmx
              hmx(nij_count) = hii(jjh)
            end if

           if ((next_nze).and.(nijcurr == jptr(jjh))) then
*              print *, nico, (hmx(ihmx),ihmx=1,nico);
              write(13) nico, (hmx(ihmx),ihmx=1,nico);
              hmx(1:nze_max) = 0.0;
              ico(1:nze_max) = 0.0;
              read_ico = .true.
           end if;

           if (next_nze) then
              nijcurr = nijcurr + 1
              nij_count = nij_count + 1
              next_nze = .false.
           end if

           if (read_ico) then
              if (jjh < ncfg) read(12) nico ,(ico(jj),jj=1,nico)
              nij_count = 1;
              jjh = jjh + 1;
              read_ico = .false.
           end if;

         end do
         ncoef = ncoef + num
         if (num .lt. ncodim) exit
      end do
      hii(1) = 0.0;
      end subroutine diag_disk_hmx


text:
!     .................................................
!             ____  _       _   ____  _____   _
!            |  _ \| |     |_| |  _ \|  ___| |_|
!            | |_) | |___   _  | |_) | |___   _
!            |  _ /|  _  | | | |  _ /|___  | | |
!            | |   | | | | | | | |    ___| | | |
!            |_|   |_| |_| |_| |_|   |_____| |_|
!     .................................................
!     PhiPsi:     a general-purpose computational
!                 mechanics program written in Fortran.
!     Website:    http://phipsi.top
!     Author:     Shi Fang from Huaiyin Institute of
!                 Technology, HuaiAn, JiangSu, China
!     Contact me: shifang@hyit.edu.cn
!     ------------------------------------------------
!     Please cite the following papers:
!     (1)Shi F, Wang X L, Liu C, Liu H, Wu H A. An
!        XFEM-based method with reduction technique
!        for modeling hydraulic fracture propagation
!        in formations containing frictional natural
!        fractures. Engineering Fracture Mechanics,
!        2017, 173: 64-90.
!     (2)Shi F, Wang X L, Liu C, Liu H, Wu H A. A
!        coupled extended finite element approach
!        for modeling hydraulic fracturing in
!        consideration of proppant. Journal of
!        Natural Gas Science and Engineering, 2016,
!        33: 885-897.
!     (3)Shi F, Wang X L, Liu C, Liu H, Wu H A. An
!        XFEM-based numerical model to calculate
!        conductivity of propped fracture considering
!        proppant transport, embedment and crushing.
!        Journal of Petroleum Science and Engineering,
!        2018, 167: 615-626..

      SUBROUTINE Vector_Unique_Int(n,n_Op_F,Vector,
     &                             Uniqued_Vec,Uniqued_n)
c     delete reduplicated element of a vector

      implicit none
      integer,intent(in):: n,n_Op_F
      integer,intent(in):: Vector(n)
      integer,intent(out)::Uniqued_Vec(n)
      integer,intent(out)::Uniqued_n

      integer i,j,k

      k = 1
      Uniqued_Vec(1) = Vector(1)

      outer: do i=2,n_Op_F
          do j=1,k
              if (Uniqued_Vec(j) .eq. Vector(i)) then
                  ! Found a match so start looking again
                  cycle outer
              end if
          end do
          ! No match found so add it to the output
          k = k + 1
          Uniqued_Vec(k) = Vector(i)
      end do outer

      Uniqued_n   = k

      return
      END SUBROUTINE Vector_Unique_Int


text:
         implicit real*8 (a-h, o-z)
         dimension a(0:500,-500:500),b(0:500,-500:500)
         dimension f(0:500,0:500),fx(0:500,0:500),fy(0:500,0:500)
         dimension fxx(0:500,0:500),fyy(0:500,0:500),fxy(0:500,0:500)
         dimension dens(0:500,0:500),dfs(0:500,0:500),dfss(0:500,0:500)

            pi=4.d0*datan(1.d0)

            open(unit=10,file='f1.dat')
            open(unit=11,file='derivatives1.dat')
            open(unit=12,file='skf1.dat')

          print 1000
 1000     format('input number of modes, number of points :',',' )
          read *, Nk,NN

          print 1001
 1001     format('input two random numbers :',',' )
          read *, xr1,yr1

c          xr1=6378836631.d0
c          yr1=77833324111.d0
          xr2=0.d0
          yr2=0.d0
          xr3=0.d0
          yr3=0.d0
c---------random number generator-----------
          ar1=0.d0
          ar2=63308.d0
          ar3=-183326.d0
          br1=86098.d0
          br2=0.d0
          br3=-539608.d0
          do 1 k1=0,Nk
             do 2 k2=-Nk,Nk
               do 3 ng=1,2
                 do 4 nrr=1,2
          xr=ar1*xr1+ar2*xr2+ar3*xr3
          yr=br1*yr1+br2*yr2+br3*yr3

          xr=dmod(xr,2147483647.d0)
          yr=dmod(yr,2145483479.d0)

          xr3=xr2
          xr2=xr1
          xr1=xr

          yr3=yr2
          yr2=yr1
          yr1=yr

          rnd=dmod((xr-yr),2147483647.d0)/2147483647.d0
          rnd=dsqrt(rnd*rnd)

          if (nrr.eq.1) then
             rnd1=rnd
               else
             rnd2=rnd
           endif
 4              continue
              grnd =dsqrt(-2.d0*dlog(rnd1))*dcos(2.d0*pi*rnd2)

            if (ng.eq.1) then
               a(k1,k2)=grnd
                else
               b(k1,k2)=grnd
            endif
 3            continue

              print *, k1,k2,a(k1,k2),b(k1,k2)

 2           continue
 1         continue
c------------------------------------------------------

           a(0,0)=0.d0
           b(0,0)=0.d0

           h=2.d0*pi/dble(NN)
            s=0.d0
           do 5 j1=0,NN-1
              do 6 j2=0,NN-1
                 x=h*dble(j1)
                 y=h*dble(j2)
                 f(j1,j2)=0.d0
                 fx(j1,j2)=0.d0
                 fy(j1,j2)=0.d0
                 fxx(j1,j2)=0.d0
                 fyy(j1,j2)=0.d0
                 fxy(j1,j2)=0.d0
                 fxxx=0.d0
                 fyyy=0.d0
                 fxxy=0.d0
                 fxyy=0.d0
                 do 7 k1=0,Nk
                    do 8 k2=-Nk,Nk
                       arg=x*dble(k1)+y*dble(k2)
                       d1=dble(k1)
                       d2=dble(k2)
               f(j1,j2)=f(j1,j2)+
     *         a(k1,k2)*dcos(arg)+b(k1,k2)*dsin(arg)
               fx(j1,j2)=fx(j1,j2)-
     *         d1*a(k1,k2)*dsin(arg)+d1*b(k1,k2)*dcos(arg)
               fy(j1,j2)=fy(j1,j2)-
     *         d2*a(k1,k2)*dsin(arg)+d2*b(k1,k2)*dcos(arg)
               fxx(j1,j2)=fxx(j1,j2)-
     *         d1*d1*a(k1,k2)*dcos(arg)-d1*d1*b(k1,k2)*dsin(arg)
               fyy(j1,j2)=fyy(j1,j2)-
     *         d2*d2*a(k1,k2)*dcos(arg)-d2*d2*b(k1,k2)*dsin(arg)
               fxy(j1,j2)=fxy(j1,j2)-
     *         d1*d2*a(k1,k2)*dcos(arg)-d1*d2*b(k1,k2)*dsin(arg)
               f111=f111+
     *         d1*d1*d1*a(k1,k2)*dsin(arg)-d1*d1*d1*b(k1,k2)*dcos(arg)
               f222=f222+
     *         d2*d2*d2*a(k1,k2)*dsin(arg)-d2*d2*d2*b(k1,k2)*dcos(arg)
               f112=f112+
     *         d1*d1*d2*a(k1,k2)*dsin(arg)-d1*d1*d2*b(k1,k2)*dcos(arg)
               f122=f122+
     *         d1*d2*d2*a(k1,k2)*dsin(arg)-d1*d2*d2*b(k1,k2)*dcos(arg)

 8                  continue
 7               continue

                 f1=fx(j1,j2)
                 f2=fy(j1,j2)
                 f11=fxx(j1,j2)
                 f22=fyy(j1,j2)
                 f12=fxy(j1,j2)

                 dfs(j1,j2)=f1*f2*(f22-f11)+(f1*f1-f2*f2)*f12

                 dfss(j1,j2)=f2*f2*(f11*f11+f1*f111+f12*f12+f2*f112)+
     *                f1*f1*(f22*f22+f2*f222+f12*f12+f1*f122)-
     *           2.d0*f1*f2*(f12*(f11+f22)+f1*f112+f2*f122)

                 s=s+f(j1,j2)*f(j1,j2)/dble(NN)/dble(NN)

                 dens(j1,j2)=f11+f22

 6            continue
 5          continue

             ss=0.d0
            do 21 j1=0,NN-1
               do 22 j2=0,NN-1
                 x=h*dble(j1)/2.d0/pi
                 y=h*dble(j2)/2.d0/pi
                  f(j1,j2)=f(j1,j2)/dsqrt(s)
                  ss=ss+f(j1,j2)*f(j1,j2)/dble(NN)/dble(NN)
                  write(10,*) x,y,f(j1,j2)
                  write(11,*) fx(j1,j2),fy(j1,j2),
     *                 fxx(j1,j2),fyy(j1,j2),fxy(j1,j2)
                  write(12,*) dens(j1,j2),dfs(j1,j2),dfss(j1,j2)
 22            continue
 21          continue

             print *, ss

             close(10)
             close(11)
             close(12)

          stop
      end


text:
\ Software and driver for the Spartan-3E Starter Kit's peripheral
\ chips connected to the SPI bus.  Contains (the beginnings of)
\ block memory words backed by the serial Flash chip.

vocabulary rekonstrukt
rekonstrukt definitions

hex
B044 constant sys-spi-base
sys-spi-base     constant sys-spi-lsb
sys-spi-base 1 + constant sys-spi-msb
sys-spi-base 2 + constant sys-spi-status
sys-spi-base 3 + constant sys-spi-config

\ serial flash routines

\ s3esk: spi flash at 30
\ 30 constant sf-address \ spi address mask
\ s3aek
01 constant sf-address \ spi address mask

: sf-initspi ( -- ) bn 0111 sys-spi-config c! ;
: sf-deselect ( -- ) 2 sf-address or sys-spi-status c! ;
: sf-waitspi ( -- )
    begin
        sys-spi-status c@
    1 and 0= until
;
: sf-xfer ( c f -- r ) \ transfer c over spi f = deselect r = result
    swap
    sys-spi-lsb c!
    if 3 else 1 then sf-address or sys-spi-status c!
    sf-waitspi
    sys-spi-lsb c@ ;

\ sf-send - send byte to flash without deselecting
: sf-send ( c -- )
    false sf-xfer drop ;

: sf-send-address ( page addr -- )
    swap sf-send
    100 /mod sf-send sf-send ;

: sf-start-reading ( page adr -- )
    03 sf-send
    sf-send-address ;

: sf-start-writing ( page adr -- )
    02 sf-send
    sf-send-address ;

: sf-read-next-word ( -- word )
    0 false sf-xfer 100 *
    0 false sf-xfer + ;

: sf-read-next-byte ( -- byte )
    0 false sf-xfer ;

: sf-rdid ( -- )
    9F sf-send
    sf-read-next-byte .
    sf-read-next-byte .
    sf-read-next-byte . cr
;

: .sf ( -- )
    05 sf-send sf-read-next-byte sf-deselect
    ." flash status: "
    dup 1 and if ." WIP " then
    dup 2 and if ." WEL " then
    dup 1C and 2 rshift ." BP:" . space
    80 and if ." SRWD " then
    cr ;

: sf-wait ( -- )
    05 sf-send
    begin
        sf-read-next-byte
    1 and 0= until
    sf-deselect ;

: sf-write-enable ( -- )
    06 sf-send sf-deselect ;

: sf-erase-page ( page -- )
    sf-write-enable
    D8 sf-send
    sf-send 0 sf-send 0 sf-send sf-deselect
    sf-wait ;

: sf-erase-bulk ( page -- )
    ." about to erase the flash, continue? "
    key 79 <> abort" aborted"
    ." erasing serial flash..."
    sf-write-enable
    C7 sf-send sf-deselect sf-wait cr ;

: sf-write-sector ( address page flash-address -- )
    sf-write-enable
    sf-start-writing
    100 0 do
        dup c@ sf-send 1+
    loop
    drop
    sf-deselect
    sf-wait ;

sf-initspi

\ Block device driver for serial flash

\ Flash can be erase in units of a page, one page is larger than one
\ block (i.e. 64K vs 1K).  Writes to the flash must be distributed to
\ the whole device, as the number of write/erase cycles per cell is
\ limited.  Thus, the number of erase cycles should be evenly
\ distributed across the device.

\ We implement a directory scheme that works as follows:

\ One page of the device is designated the directory page.  For every
\ allocatable physical block, there is one two byte entry in the
\ directory page that indicates the logical block number to which
\ the physical block is allocated.  FFFF signifies an unused block.
\ To find the physical block number for a given logical block number,
\ the directory is sequentially scanned until an entry with the
\ required block number is found or FFFF is encountered.  When
\ FFFF has been encountered, the logical block has not previously
\ been allocated.
\ In order to write a logical block, the first unused (FFFF) entry
\ in the directory is overwritten with the logical block number.

\ Once the directory is completely filled, unused entries need to be
\ reclaimed.  This is done by copying all entries from the currently
\ active directory page to the secondary directory page and than erasing
\ the currently active directory page.

\ When all blocks have been allocated, a compaction cycle must be run.

\ When starting, the primary and secondary directory page are scanned.
\ If both pages are empty, the primary page becomes the active page.
\ If one of them is empty, the non-empty page becomes the active page.
\ It is an error if none of the directory pages are empty.

\ Serial flash definitions, 16 MBit flash

decimal
32 constant pages     \ Number of pages
64 constant page-size \ Size of one page in 1K blocks

\ Number of data pages
30 constant data-page-count
\ Page addresses of primary and secondary directory page
31 constant primary-directory
30 constant secondary-directory
hex

variable active-directory-page

: page-empty? ( n -- f )
    03 sf-send
    sf-send \ send page number
    0 sf-send 0 sf-send
    true
    FFFF 0 do
        0 false sf-xfer
        FF = 0= if
            drop false leave
        then
    loop
    0 true sf-xfer drop ;

: ?internal-error ( f -- )
    true abort" Internal error, both directory pages occupied" ;

: find-active-directory-page ( -- )
    secondary-directory page-empty? if
        primary-directory active-directory-page !
    else
        primary-directory page-empty? 0= ?internal-error
        secondary-directory active-directory-page !
    then ;

: sf-directory-addr ( phys-block-no -- page addr )
    active-directory-page @ swap 1 lshift ;

: sf-open-directory ( -- )
    0 sf-directory-addr sf-start-reading ;

: scanning-directory ( -- end begin )
    data-page-count page-size *
    0 ;

\ Find physical block number - n is the logical block number to
\ find, returns -1 if not found
: find-block ( log-block-no -- phys-block-no )
    sf-open-directory
    FFFF
    scanning-directory do
        sf-read-next-word
        2 pick = if
            drop i
        then
    loop
    sf-deselect
    swap drop ;

: ?no-free-block ( f -- )
    abort" no more free blocks" ;

\ Find free physical block
: find-free-block ( -- phys-block-no )
    sf-open-directory
    FFFF
    scanning-directory do
        sf-read-next-word
        FFFF = if
            drop i leave
        then
    loop
    sf-deselect
    dup FFFF = ?no-free-block ;

\ Show usage statistics
: .sf-stat ( -- )
    sf-open-directory
    0 0
    scanning-directory do
        sf-read-next-word
        FFFF = if
            1+
        else
            swap 1+ swap
        then
    loop
    sf-deselect
    ." free: " . ." in use: " . cr ;

: ?entry-in-use ( f -- )
    abort" can't overwrite directory entry that is in use" ;

: assign-block ( log-block-no phys-block-no -- )
    sf-write-enable
    sf-directory-addr sf-start-writing
    100 /mod sf-send sf-send sf-deselect sf-wait ;

variable updated
variable scr
variable block-buf 400 allot

: .block ( -- )
    block-buf 100 dump \ only print first 256 bytes
    ;

: sf-read-block ( page addr -- )
    sf-start-reading
    block-buf 400 + block-buf do
        sf-read-next-byte i c!
    loop
    sf-deselect ;

: block-to-flash ( physblock -- page addr )
    2 lshift 100 /mod swap 8 lshift ;

: block ( block -- addr )
    find-block
    dup FFFF = if
        drop
        block-buf 400 0 fill
    else
        block-to-flash sf-read-block
    then
    block-buf ;

: update ( -- )
    1 updated ! ;

: flush ( -- )
    updated @ if
        find-free-block
        400 0 do
            block-buf i +
            over block-to-flash i +
            sf-write-sector
        100 +loop
        scr @ swap ." assigning block " .s assign-block
        false updated !
    then ;

\ ADC tests

hex
: adc-initspi ( -- ) bn 1111 sys-spi-config c! ;
: adc-read ( -- word )
    23 sys-spi-status c!
    sys-spi-msb @ 8 rshift
    sys-spi-lsb @ or ;
: adc-loop ( -- ) begin adc-read . cr key 20 <> until ;

max_stars_repo_name: hanshuebner/rekonstrukt

max_issues_repo_name: hanshuebner/rekonstrukt

max_forks_repo_name: hanshuebner/rekonstrukt


text:
Healing Hands Massage is the Massage Services massage practice of owner Users/KarenKhamashta Karen S. Khamashta, CMT.  Operating since 1982, it is a well established local business.

Karens expertise is her wellstudied understanding of muscle physiology and anatomical structure that enables her to effectively help those suffering from acute or chronic tissue pain, drawing upon muscle release bodywork techniques and restorative medical massage and 23 years of experience.

  From therapy to pure relaxation, Davis offers a variety of Massage Services

20120317 15:31:36 nbsp I have been going to Karen Khamashta for regular massage therapy since 2008 when my previous massage therapist discontinued her practice.  Karen specializes in trigger point massage which involves painlessly releasing muscle knots that are believed to be a source of pain experienced elsewhere in the body. It took me a while get accustomed to her style as I was used to a no interaction Swedish style relaxation approach, and wasn’t looking specifically for pain relief. However, after ten minutes on the table I now feel the tension begin to melt away. By the end of the massage I am completely relaxed and any incipient pain in my body is gone. Karen’s office is located in a quiet neighborhood in Davis.  She uses a heated pad on the table so I am never cold, even in the middle of winter  something I really appreciate. And her fee  $85 for an hour and a half – is very reasonable. Whether you are looking for relief from chronic pain or just relief from everyday stress, Karen can help you.

 Users/MaryFontaine

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


!***********************************/
!	Name:ローカル格子に対応した補助構造格子の生成
!	Alias:OMakeAuxiliaryStructuredGrid
!	Description:ローカル格子移動計算を効率化するための補助格子を生成する
!	Type:
!	Input:UG,MG
!	Output:
!	Note:補助構造格子の各セル内にセル中心が存在するローカル格子を登録し，相互に検索できるようにする
!	Author:Akitaka Toyota
!	Date:2017.11.25
!	Update:2017.12.29
!	Other:適切なUG,MGさえ入れれば補助格子が返るため，複数格子用の変更は不要
!***********************************/
subroutine OMakeAuxiliaryStructuredGrid(UG,RC,ASG)
    use StructVar_Mod
    use LoopVar_Mod
    implicit none
    type(UnstructuredGrid), intent(in) :: UG
    type(RelativeCoordinate), intent(in) :: RC
    type(AuxiliaryStructuredGrid), intent(inout) :: ASG

    double precision :: AverageVolume, TmpWidth
!    type(Structured2Unstructured), pointer :: StU
    write(6,*) "Make Auxiliary Grid for Overset"
    call MakeAuxiliaryGrid
    call CorrespondGrid
    !call CheckAuxiliaryGrid

    !do iCell=1, ASG%TotalCell

    !StU => ASG%AS2U(iCell)%S2U
    !if (associated(StU)) write(6,*) "Target",iCell
    !do while(associated(StU))
    !    write(6,*) StU%Cell
    !    StU => StU%Next
    !end do
    !end do
    !stop

return
contains

    subroutine MakeAuxiliaryGrid
    implicit none
    double precision :: SmallIncrement
    !for 2D-Grid

        SmallIncrement = 10.0d0**(-9)
    !この値は本当に適当で，0より大きい十分に小さな値であれば正直なんでもいい
    !この値の幅だけローカル格子からはみ出す無駄な領域が増えるため，計算中に消えない程度にできるだけ小さく取ればそれで十分

        ASG%Bound(1,1) = minval(RC%Cell(:,1),1) - SmallIncrement !まったく同じ座標にしてしまうとグローバル番号を検索する式が不正な値を取り，配列がオーバーフローしてしまう
        ASG%Bound(2,1) = maxval(RC%Cell(:,1),1) + SmallIncrement
        ASG%Bound(1,2) = minval(RC%Cell(:,2),1) - SmallIncrement
        ASG%Bound(2,2) = maxval(RC%Cell(:,2),1) + SmallIncrement
        ASG%Bound(1,3) = minval(RC%Cell(:,3),1) - SmallIncrement
        ASG%Bound(2,3) = maxval(RC%Cell(:,3),1) + SmallIncrement
        !width fix to integer
        AverageVolume = sum(UG%GM%Volume(:)) / dble(UG%GI%RealCells) !要素の平均体積を計算する
        TmpWidth = 2.0d0* sqrt(AverageVolume)

        ASG%CellNumber(1:2) = int((ASG%Bound(2,1:2) - ASG%Bound(1,1:2))/TmpWidth)
        ASG%CellNumber(3) = 1

        ASG%Width(:) = (ASG%Bound(2,:) - ASG%Bound(1,:))/dble(ASG%CellNumber(:))
        ASG%TotalCell = ASG%CellNumber(1) * ASG%CellNumber(2) * ASG%CellNumber(3)

        return
    end subroutine MakeAuxiliaryGrid

    subroutine CorrespondGrid
    use FrequentOperation
    implicit none
    !type(RelativeCoordinate), intent(in) :: RC
    !type(AuxiliaryStructuredGrid), intent(inout) :: ASG

    type(Structured2Unstructured), pointer :: Temp
    integer :: iGlobalCellNumber

    allocate(ASG%RelatedU2AS(UG%GI%AllCells,3)) !(非構造要素番号,xyz)で補助格子の座標軸に沿った番号を返す !この変数いらない説が濃厚 !Local変数として破棄して大丈夫なやーつっぽい
    allocate(ASG%ContentsAS2U(ASG%TotalCell)) !(補助格子のグローバル番号)で補助格子が保有する非構造要素数が返る
    allocate(ASG%AS2U(ASG%TotalCell)) !(補助格子のグローバル番号)%S2U%Cellで内包セル，%S2U%Nextで次の内包セル

        ASG%ContentsAS2U = 0

        do iCell=1, UG%GI%RealCells !すべての非構造実セルについて
            !Auxiliary格子の座標は配列の添字と対応しているため，
            !(非構造格子の座標r,Auxiliary格子の座標r0，構造格子負端からのセル数n，格子幅dxについて
            !r-r0 > n*dxなる関係があることから，nについて解くと以下のような式になる(※2次元用で書いたためzは常に固定値になる)

            call GetLocalNumber2D(RC%Cell(iCell,1:2),ASG%Width(1:2),ASG%Bound(1,1:2),ASG%RelatedU2AS(iCell,1:2))
            ASG%RelatedU2AS(iCell,3) = 1

            !取扱いの利便性を上げるためグローバルセル番号(詳細は↓)に変換する．
            !(i,j,k) <=>!2D:(j-1)*(Lx)+i !3D:(k-1)*(Ly*Lx)+(j-1)*(jmax-1)+i
            iGlobalCellNumber = Local2Global(ASG%RelatedU2AS(iCell,1),ASG%RelatedU2AS(iCell,2),ASG%CellNumber(1))

            !Auxiliary格子ijkが非構造格子から参照された回数を記録
            ASG%ContentsAS2U(iGlobalCellNumber) = ASG%ContentsAS2U(iGlobalCellNumber) + 1

            if(ASG%ContentsAS2U(iGlobalCellNumber) == 1) then
                nullify(ASG%AS2U(iGlobalCellNumber)%S2U) !リストの初期化
            end if
            !最後に構造格子に所属する非構造格子への対応を記録する
            allocate(Temp) !新規ノードの確保
            Temp%Cell = iCell !ノードへの記録
            Temp%Next => ASG%AS2U(iGlobalCellNumber)%S2U !前のノードへのリンク作成
            ASG%AS2U(iGlobalCellNumber)%S2U => Temp !データの保存

        end do

        return
    end subroutine CorrespondGrid

    subroutine CheckAuxiliaryGrid
    use LoopVar_Mod
    implicit none

    open(unit = 1, file ="Auxiliary1.vtk", status = 'unknown')
        write(1,"('# vtk DataFile Version 3.0')")
        write(1,*) "AuxiliaryTest"
        write(1,"('ASCII')")
        write(1,"('DATASET UNSTRUCTURED_GRID')")
        write(1,"('POINTS ',(1x,i7),' double')") UG%GI%Points

        do iPoint=1, UG%GI%Points
            !write(1,"(3(1x,f22.17))") (UG%CD%Point(iPoint,iLoop),iLoop=1,3) !Adjust Point Number
            write(1,"(3(1x,f22.17))") UG%CD%Point(iPoint,1),UG%CD%Point(iPoint,2),UG%CD%Point(iPoint,3)+1 !for OverSet
        end do
        write(1,*) ""

        write(1,"('CELLS ',2(1x,i7))") UG%GI%RealCells,UG%GI%RealCells*4
        do iCell=1, UG%GI%RealCells
                write(1,"(4(1x,i7))") 3,(UG%Tri%Point(iCell,iLoop)-1,iLoop=1,3) !Adjusted Point Number
        end do
        write(1,*) ""

        write(1,"('CELL_TYPES ',(1x,i7))") UG%GI%RealCells
        do iCell=1, UG%GI%RealCells
            write(1,"(1x,i1)") 5
        end do

        write(1,"('CELL_DATA ',i7)") UG%GI%RealCells

        write(1,"('SCALARS AuxiliaryNumber float')")
        write(1,"('LOOKUP_TABLE default')")
        do iCell=1, UG%GI%RealCells
            write(1,"((2x,i6))") (ASG%RelatedU2AS(iCell,2)-1)*ASG%CellNumber(1) + ASG%RelatedU2AS(iCell,1)
            !write(6,*) (ASG%RelatedU2AS(iCell,2)-1)*ASG%CellNumber(1) + ASG%RelatedU2AS(iCell,1)
        end do
    close(1)

    open(unit=1,file="Auxiliary2.vtk",status='unknown')
        write(1,"('# vtk DataFile Version 3.0')")
        write(1,"('Auxiliary2')")
        write(1,"('ASCII')")
        write(1,"('DATASET RECTILINEAR_GRID')")
        write(1,"('DIMENSIONS ',3(1x,i5))") ASG%CellNumber(1)+1,ASG%CellNumber(2)+1,1

!座標情報
        write(1,"('X_COORDINATES ',i5,' float')") ASG%CellNumber(1)+1
            do iCenterX=1, ASG%CellNumber(1)+1
                write(1,"(f22.17)") ASG%Bound(1,1) + ASG%Width(1) * (dble(iCenterX)-1.0d0)
            end do

        write(1,"('Y_COORDINATES ',i5,' float')") ASG%CellNumber(2)+1
        do iCenterY=1, ASG%CellNumber(2)+1
            write(1,"(f22.17)") ASG%Bound(1,2) + ASG%Width(2) * (dble(iCenterY)-1.0d0)
        end do

        write(1,"('Z_COORDINATES ',i5,' float')") 1
        write(1,"(f22.17)") 0.0

        write(1,"('CELL_DATA ',i7)") (ASG%CellNumber(1)) * (ASG%CellNumber(2))

        write(1,"('SCALARS Contents float')")
        write(1,"('LOOKUP_TABLE default')")
            do iCenterY=1,ASG%CellNumber(2)
                do iCenterX=1,ASG%CellNumber(1)
                    write(1,"((2x,i6))") ASG%ContentsAS2U((iCenterY-1)*ASG%CellNumber(1) + iCenterX)
                end do
            end do

        write(1,"('SCALARS GlobalNumber float')")
        write(1,"('LOOKUP_TABLE default')")
            do iCenterY=1,ASG%CellNumber(2)
                do iCenterX=1,ASG%CellNumber(1)
                    write(1,"((2x,i6))") (iCenterY-1)*ASG%CellNumber(1) + iCenterX
                end do
            end do

        close(1)

    return
    end subroutine CheckAuxiliaryGrid

end subroutine OMakeAuxiliaryStructuredGrid


text:
module file_with_ext_sub_test
  use fruit
contains
  subroutine test_uppercase
    use abc_module

    character(len = 1) :: ch_result
    call uppercase("a", ch_result)
    call assert_equals("A", ch_result)

    call uppercase("#", ch_result)
    call assert_equals("#", ch_result)

  end subroutine test_uppercase
end module file_with_ext_sub_test

max_stars_repo_name: thomasms/fortran-examples

max_issues_repo_name: thomasms/fortran-examples

max_forks_repo_name: thomasms/fortran-examples


text:
 SUBROUTINE POLATES3(IPOPT,IGDTNUMI,IGDTMPLI,IGDTLENI, &
                     IGDTNUMO,IGDTMPLO,IGDTLENO, &
                     MI,MO,KM,IBI,LI,GI, &
                     NO,RLAT,RLON,IBO,LO,GO,IRET)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!
! $Revision: 74685 $
!
! SUBPROGRAM:  POLATES3   INTERPOLATE SCALAR FIELDS (BUDGET)
!   PRGMMR: IREDELL       ORG: W/NMC23       DATE: 96-04-10
!
! ABSTRACT: THIS SUBPROGRAM PERFORMS BUDGET INTERPOLATION
!           FROM ANY GRID TO ANY GRID (OR TO RANDOM STATION
!           POINTS) FOR SCALAR FIELDS.  THE ALGORITHM
!           SIMPLY COMPUTES (WEIGHTED) AVERAGES OF
!           BILINEARLY INTERPOLATED POINTS ARRANGED IN A SQUARE BOX
!           CENTERED AROUND EACH OUTPUT GRID POINT AND STRETCHING
!           NEARLY HALFWAY TO EACH OF THE NEIGHBORING GRID POINTS.
!           OPTIONS ALLOW CHOICES OF NUMBER OF POINTS IN EACH RADIUS
!           FROM THE CENTER POINT (IPOPT(1)) WHICH DEFAULTS TO 2
!           (IF IPOPT(1)=-1) MEANING THAT 25 POINTS WILL BE AVERAGED;
!           FURTHER OPTIONS ARE THE RESPECTIVE WEIGHTS FOR THE RADIUS
!           POINTS STARTING AT THE CENTER POINT (IPOPT(2:2+IPOPT(1))
!           WHICH DEFAULTS TO ALL 1 (IF IPOPT(1)=-1 OR IPOPT(2)=-1).
!           A SPECIAL INTERPOLATION IS DONE IF IPOPT(2)=-2.
!           IN THIS CASE, THE BOXES STRETCH NEARLY ALL THE WAY TO
!           EACH OF THE NEIGHBORING GRID POINTS AND THE WEIGHTS
!           ARE THE ADJOINT OF THE BILINEAR INTERPOLATION WEIGHTS.
!           THIS CASE GIVES QUASI-SECOND-ORDER BUDGET INTERPOLATION.
!           ANOTHER OPTION IS THE MINIMUM PERCENTAGE FOR MASK,
!           I.E. PERCENT VALID INPUT DATA REQUIRED TO MAKE OUTPUT DATA,
!           (IPOPT(3+IPOPT(1)) WHICH DEFAULTS TO 50 (IF -1).
!           IN CASES WHERE THERE IS NO OR INSUFFICIENT VALID INPUT DATA,
!           THE USER MAY CHOOSE TO SEARCH FOR THE NEAREST VALID DATA.
!           THIS IS INVOKED BY SETTING IPOPT(20) TO THE WIDTH OF
!           THE SEARCH SQUARE. THE DEFAULT IS 1 (NO SEARCH).  SQUARES ARE
!           SEARCHED FOR VALID DATA IN A SPIRAL PATTERN
!           STARTING FROM THE CENTER.  NO SEARCHING IS DONE WHERE
!           THE OUTPUT GRID IS OUTSIDE THE INPUT GRID.
!           ONLY HORIZONTAL INTERPOLATION IS PERFORMED.
!
!           THE CODE RECOGNIZES THE FOLLOWING PROJECTIONS, WHERE
!           "IGDTNUMI/O" IS THE GRIB 2 GRID DEFINTION TEMPLATE NUMBER
!           FOR THE INPUT AND OUTPUT GRIDS, RESPECTIVELY:
!             (IGDTNUMI/O=00) EQUIDISTANT CYLINDRICAL
!             (IGDTNUMI/O=01) ROTATED EQUIDISTANT CYLINDRICAL. "E" AND
!                             NON-"E" STAGGERED
!             (IGDTNUMI/O=10) MERCATOR CYLINDRICAL
!             (IGDTNUMI/O=20) POLAR STEREOGRAPHIC AZIMUTHAL
!             (IGDTNUMI/O=30) LAMBERT CONFORMAL CONICAL
!             (IGDTNUMI/O=40) GAUSSIAN CYLINDRICAL
!           AS AN ADDED BONUS THE NUMBER OF OUTPUT GRID POINTS
!           AND THEIR LATITUDES AND LONGITUDES ARE ALSO RETURNED.
!           ON THE OTHER HAND, THE OUTPUT CAN BE A SET OF STATION POINTS
!           IF IGDTNUMO=IGDTNUMO-255, IN WHICH CASE THE NUMBER OF POINTS
!           AND THEIR LATITUDES AND LONGITUDES MUST BE INPUT.
!
!           INPUT BITMAPS WILL BE INTERPOLATED TO OUTPUT BITMAPS.
!           OUTPUT BITMAPS WILL ALSO BE CREATED WHEN THE OUTPUT GRID
!           EXTENDS OUTSIDE OF THE DOMAIN OF THE INPUT GRID.
!           THE OUTPUT FIELD IS SET TO 0 WHERE THE OUTPUT BITMAP IS OFF.
!
! PROGRAM HISTORY LOG:
!   96-04-10  IREDELL
! 1999-04-08  IREDELL  SPLIT IJKGDS INTO TWO PIECES
! 1999-04-08  IREDELL  ADDED BILINEAR OPTION IPOPT(2)=-2
! 2001-06-18  IREDELL  INCLUDE MINIMUM MASK PERCENTAGE OPTION
! 2006-01-04  GAYNO    ADDED OPTION TO DO SUBSECTION OF OUTPUT GRID.
!                      ADDED SPIRAL SEARCH OPTION.
! 2015-01-27  GAYNO    REPLACE CALLS TO GDSWIZ WITH NEW MERGED
!                      VERSION OF GDSWZD.
! 2015-07-13  GAYNO    CONVERT TO GRIB 2. REPLACE GRIB 1 KGDS ARRAYS
!                      WITH GRIB 2 GRID DEFINITION TEMPLATE ARRAYS.
!
! USAGE:    CALL POLATES3(IPOPT,IGDTNUMI,IGDTMPLI,IGDTLENI, &
!                         IGDTNUMO,IGDTMPLO,IGDTLENO, &
!                         MI,MO,KM,IBI,LI,GI, &
!                         NO,RLAT,RLON,IBO,LO,GO,IRET)
!
!   INPUT ARGUMENT LIST:
!     IPOPT    - INTEGER (20) INTERPOLATION OPTIONS
!                IPOPT(1) IS NUMBER OF RADIUS POINTS
!                (DEFAULTS TO 2 IF IPOPT(1)=-1);
!                IPOPT(2:2+IPOPT(1)) ARE RESPECTIVE WEIGHTS
!                (DEFAULTS TO ALL 1 IF IPOPT(1)=-1 OR IPOPT(2)=-1).
!                IPOPT(3+IPOPT(1)) IS MINIMUM PERCENTAGE FOR MASK
!                (DEFAULTS TO 50 IF IPOPT(3+IPOPT(1)=-1)
!     IGDTNUMI - INTEGER GRID DEFINITION TEMPLATE NUMBER - INPUT GRID.
!                CORRESPONDS TO THE GFLD%IGDTNUM COMPONENT OF THE
!                NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE:
!                  00 - EQUIDISTANT CYLINDRICAL
!                  01 - ROTATED EQUIDISTANT CYLINDRICAL.  "E"
!                       AND NON-"E" STAGGERED
!                  10 - MERCATOR CYCLINDRICAL
!                  20 - POLAR STEREOGRAPHIC AZIMUTHAL
!                  30 - LAMBERT CONFORMAL CONICAL
!                  40 - GAUSSIAN EQUIDISTANT CYCLINDRICAL
!     IGDTMPLI - INTEGER (IGDTLENI) GRID DEFINITION TEMPLATE ARRAY -
!                INPUT GRID. CORRESPONDS TO THE GFLD%IGDTMPL COMPONENT
!                OF THE NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
!                (SECTION 3 INFO).  SEE COMMENTS IN ROUTINE
!                IPOLATES FOR COMPLETE DEFINITION.
!     IGDTLENI - INTEGER NUMBER OF ELEMENTS OF THE GRID DEFINITION
!                TEMPLATE ARRAY - INPUT GRID.  CORRESPONDS TO THE GFLD%IGDTLEN
!                COMPONENT OF THE NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
!     IGDTNUMO - INTEGER GRID DEFINITION TEMPLATE NUMBER - OUTPUT GRID.
!                CORRESPONDS TO THE GFLD%IGDTNUM COMPONENT OF THE
!                NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE. IGDTNUMO=IGDTNUM-255
!                MEANS INTERPOLATE TO RANDOM STATION POINTS.
!                OTHERWISE, SAME DEFINITION AS "IGDTNUMI".
!     IGDTMPLO - INTEGER (IGDTLENO) GRID DEFINITION TEMPLATE ARRAY -
!                OUTPUT GRID. CORRESPONDS TO THE GFLD%IGDTMPL COMPONENT
!                OF THE NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
!                (SECTION 3 INFO).  SEE COMMENTS IN ROUTINE
!                IPOLATES FOR COMPLETE DEFINITION.
!     IGDTLENO - INTEGER NUMBER OF ELEMENTS OF THE GRID DEFINITION
!                TEMPLATE ARRAY - OUTPUT GRID.  CORRESPONDS TO THE GFLD%IGDTLEN
!                COMPONENT OF THE NCEP G2 LIBRARY GRIDMOD DATA STRUCTURE.
!     MI       - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1
!                OR DIMENSION OF INPUT GRID FIELDS IF KM=1
!     MO       - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1
!                OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1
!     KM       - INTEGER NUMBER OF FIELDS TO INTERPOLATE
!     IBI      - INTEGER (KM) INPUT BITMAP FLAGS
!     LI       - LOGICAL*1 (MI,KM) INPUT BITMAPS (IF SOME IBI(K)=1)
!     GI       - REAL (MI,KM) INPUT FIELDS TO INTERPOLATE
!     RLAT     - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF IGDTNUMO=IGDTNUMO-255)
!     RLON     - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF IGDTNUMO=IGDTNUMO-255)
!
!   OUTPUT ARGUMENT LIST:
!     NO       - INTEGER NUMBER OF OUTPUT POINTS
!     RLAT     - REAL (MO) OUTPUT LATITUDES IN DEGREES (ONLY IF IGDTNUMO>=0)
!     RLON     - REAL (MO) OUTPUT LONGITUDES IN DEGREES (ONLY IF IGDTNUMO>=0)
!     IBO      - INTEGER (KM) OUTPUT BITMAP FLAGS
!     LO       - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT)
!     GO       - REAL (MO,KM) OUTPUT FIELDS INTERPOLATED
!     IRET     - INTEGER RETURN CODE
!                0    SUCCESSFUL INTERPOLATION
!                2    UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP
!                3    UNRECOGNIZED OUTPUT GRID
!                32   INVALID BUDGET METHOD PARAMETERS
!
! SUBPROGRAMS CALLED:
!   GDSWZD       GRID DESCRIPTION SECTION WIZARD
!   IJKGDS0      SET UP PARAMETERS FOR IJKGDS1
!   IJKGDS1      RETURN FIELD POSITION FOR A GIVEN GRID POINT
!   POLFIXS      MAKE MULTIPLE POLE SCALAR VALUES CONSISTENT
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!
!$$$
 USE GDSWZD_MOD

 IMPLICIT NONE
!
 INTEGER,    INTENT(IN   )     :: IGDTNUMI, IGDTLENI
 INTEGER,    INTENT(IN   )     :: IGDTMPLI(IGDTLENI)
 INTEGER,    INTENT(IN   )     :: IGDTNUMO, IGDTLENO
 INTEGER,    INTENT(IN   )     :: IGDTMPLO(IGDTLENO)
 INTEGER,    INTENT(IN   )     :: IBI(KM), IPOPT(20)
 INTEGER,    INTENT(IN   )     :: KM, MI, MO
 INTEGER,    INTENT(  OUT)     :: IBO(KM), IRET, NO
!
 LOGICAL*1,  INTENT(IN   )     :: LI(MI,KM)
 LOGICAL*1,  INTENT(  OUT)     :: LO(MO,KM)
!
 REAL,       INTENT(IN   )     :: GI(MI,KM)
 REAL,       INTENT(INOUT)     :: RLAT(MO), RLON(MO)
 REAL,       INTENT(  OUT)     :: GO(MO,KM)
!
 REAL,       PARAMETER         :: FILL=-9999.
!
 INTEGER                       :: IJKGDS1, I1, J1, I2, J2, IB, JB
 INTEGER                       :: IJKGDSA(20), IX, JX, IXS, JXS
 INTEGER                       :: K, KXS, KXT, IGDTNUMO2
 INTEGER                       :: LB, LSW, MP, MSPIRAL, MX
 INTEGER                       :: N, NB, NB1, NB2, NB3, NB4, NV, NX
 INTEGER                       :: N11(MO),N21(MO),N12(MO),N22(MO)
!
 REAL                          :: GB, LAT(1), LON(1)
 REAL                          :: PMP, RB2, RLOB(MO), RLAB(MO), WB
 REAL                          :: W11(MO), W21(MO), W12(MO), W22(MO)
 REAL                          :: WO(MO,KM), XF, YF, XI, YI, XX, YY
 REAL                          :: XPTS(MO),YPTS(MO),XPTB(MO),YPTB(MO)
 REAL                          :: XXX(1), YYY(1)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES.
!  DO SUBSECTION OF GRID IF KGDSO(1) IS SUBTRACTED FROM 255.
 IRET=0
 IF(IGDTNUMO.GE.0) THEN
   CALL GDSWZD(IGDTNUMO,IGDTMPLO,IGDTLENO, 0,MO,FILL,XPTS,YPTS,RLON,RLAT,NO)
   IF(NO.EQ.0) IRET=3
 ELSE
   IGDTNUMO2=255+IGDTNUMO
   CALL GDSWZD(IGDTNUMO2,IGDTMPLO,IGDTLENO,-1,MO,FILL,XPTS,YPTS,RLON,RLAT,NO)
   IF(NO.EQ.0) IRET=3
 ENDIF
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  SET PARAMETERS
 IF(IPOPT(1).GT.16) IRET=32
 MSPIRAL=MAX(IPOPT(20),1)
 NB1=IPOPT(1)
 IF(NB1.EQ.-1) NB1=2
 IF(IRET.EQ.0.AND.NB1.LT.0) IRET=32
 LSW=1
 IF(IPOPT(2).EQ.-2) LSW=2
 IF(IPOPT(1).EQ.-1.OR.IPOPT(2).EQ.-1) LSW=0
 IF(IRET.EQ.0.AND.LSW.EQ.1.AND.NB1.GT.15) IRET=32
 MP=IPOPT(3+IPOPT(1))
 IF(MP.EQ.-1.OR.MP.EQ.0) MP=50
 IF(MP.LT.0.OR.MP.GT.100) IRET=32
 PMP=MP*0.01
 IF(IRET.EQ.0) THEN
   NB2=2*NB1+1
   RB2=1./NB2
   NB3=NB2*NB2
   NB4=NB3
   IF(LSW.EQ.2) THEN
     RB2=1./(NB1+1)
     NB4=(NB1+1)**4
   ELSEIF(LSW.EQ.1) THEN
     NB4=IPOPT(2)
     DO IB=1,NB1
       NB4=NB4+8*IB*IPOPT(2+IB)
     ENDDO
   ENDIF
 ELSE
   NB3=0
   NB4=1
 ENDIF
 DO K=1,KM
   DO N=1,NO
     GO(N,K)=0.
     WO(N,K)=0.
   ENDDO
 ENDDO
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  LOOP OVER SAMPLE POINTS IN OUTPUT GRID BOX
 CALL IJKGDS0(IGDTNUMI,IGDTMPLI,IGDTLENI,IJKGDSA)
 DO NB=1,NB3
!  LOCATE INPUT POINTS AND COMPUTE THEIR WEIGHTS
   JB=(NB-1)/NB2-NB1
   IB=NB-(JB+NB1)*NB2-NB1-1
   LB=MAX(ABS(IB),ABS(JB))
   WB=1
   IF(LSW.EQ.2) THEN
     WB=(NB1+1-ABS(IB))*(NB1+1-ABS(JB))
   ELSEIF(LSW.EQ.1) THEN
     WB=IPOPT(2+LB)
   ENDIF
   IF(WB.NE.0) THEN
!$OMP PARALLEL DO PRIVATE(N) SCHEDULE(STATIC)
     DO N=1,NO
       XPTB(N)=XPTS(N)+IB*RB2
       YPTB(N)=YPTS(N)+JB*RB2
     ENDDO
!$OMP END PARALLEL DO
     CALL GDSWZD(IGDTNUMO,IGDTMPLO,IGDTLENO, 1,NO,FILL,XPTB,YPTB,RLOB,RLAB,NV)
     CALL GDSWZD(IGDTNUMI,IGDTMPLI,IGDTLENI,-1,NO,FILL,XPTB,YPTB,RLOB,RLAB,NV)
     IF(IRET.EQ.0.AND.NV.EQ.0.AND.LB.EQ.0) IRET=2
!$OMP PARALLEL DO PRIVATE(N,XI,YI,I1,I2,J1,J2,XF,YF) SCHEDULE(STATIC)
     DO N=1,NO
       XI=XPTB(N)
       YI=YPTB(N)
       IF(XI.NE.FILL.AND.YI.NE.FILL) THEN
         I1=XI
         I2=I1+1
         J1=YI
         J2=J1+1
         XF=XI-I1
         YF=YI-J1
         N11(N)=IJKGDS1(I1,J1,IJKGDSA)
         N21(N)=IJKGDS1(I2,J1,IJKGDSA)
         N12(N)=IJKGDS1(I1,J2,IJKGDSA)
         N22(N)=IJKGDS1(I2,J2,IJKGDSA)
         IF(MIN(N11(N),N21(N),N12(N),N22(N)).GT.0) THEN
           W11(N)=(1-XF)*(1-YF)
           W21(N)=XF*(1-YF)
           W12(N)=(1-XF)*YF
           W22(N)=XF*YF
         ELSE
           N11(N)=0
           N21(N)=0
           N12(N)=0
           N22(N)=0
         ENDIF
       ELSE
         N11(N)=0
         N21(N)=0
         N12(N)=0
         N22(N)=0
       ENDIF
     ENDDO
!$OMP END PARALLEL DO
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  INTERPOLATE WITH OR WITHOUT BITMAPS
!$OMP PARALLEL DO PRIVATE(K,N,GB) SCHEDULE(STATIC)
     DO K=1,KM
       DO N=1,NO
         IF(N11(N).GT.0) THEN
           IF(IBI(K).EQ.0) THEN
             GB=W11(N)*GI(N11(N),K)+W21(N)*GI(N21(N),K) &
                +W12(N)*GI(N12(N),K)+W22(N)*GI(N22(N),K)
             GO(N,K)=GO(N,K)+WB*GB
             WO(N,K)=WO(N,K)+WB
           ELSE
             IF(LI(N11(N),K)) THEN
               GO(N,K)=GO(N,K)+WB*W11(N)*GI(N11(N),K)
               WO(N,K)=WO(N,K)+WB*W11(N)
             ENDIF
             IF(LI(N21(N),K)) THEN
               GO(N,K)=GO(N,K)+WB*W21(N)*GI(N21(N),K)
               WO(N,K)=WO(N,K)+WB*W21(N)
             ENDIF
             IF(LI(N12(N),K)) THEN
               GO(N,K)=GO(N,K)+WB*W12(N)*GI(N12(N),K)
               WO(N,K)=WO(N,K)+WB*W12(N)
             ENDIF
             IF(LI(N22(N),K)) THEN
               GO(N,K)=GO(N,K)+WB*W22(N)*GI(N22(N),K)
               WO(N,K)=WO(N,K)+WB*W22(N)
             ENDIF
           ENDIF
         ENDIF
       ENDDO
     ENDDO
!$OMP END PARALLEL DO
   ENDIF
 ENDDO   ! sub-grid points
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  COMPUTE OUTPUT BITMAPS AND FIELDS
! KM is often 1 .. do not do OMP PARALLEL DO here
 KM_LOOP : DO K=1,KM
   IBO(K)=IBI(K)
!$OMP PARALLEL DO PRIVATE(N,LAT,LON,XXX,YYY,NV,XX,YY,IXS,JXS,MX,KXS,KXT,IX,JX,NX) SCHEDULE(STATIC)
   N_LOOP : DO N=1,NO
     LO(N,K)=WO(N,K).GE.PMP*NB4
     IF(LO(N,K)) THEN
       GO(N,K)=GO(N,K)/WO(N,K)
     ELSEIF (MSPIRAL.GT.1) THEN
       LAT(1)=RLAT(N)
       LON(1)=RLON(N)
       CALL GDSWZD(IGDTNUMI,IGDTMPLI,IGDTLENI,-1,1,FILL,XXX,YYY,LON,LAT,NV)
       XX=XXX(1)
       YY=YYY(1)
       IF(NV.EQ.1)THEN
         I1=NINT(XX)
         J1=NINT(YY)
         IXS=SIGN(1.,XX-I1)
         JXS=SIGN(1.,YY-J1)
         SPIRAL_LOOP : DO MX=2,MSPIRAL**2
           KXS=SQRT(4*MX-2.5)
           KXT=MX-(KXS**2/4+1)
           SELECT CASE(MOD(KXS,4))
           CASE(1)
             IX=I1-IXS*(KXS/4-KXT)
             JX=J1-JXS*KXS/4
           CASE(2)
             IX=I1+IXS*(1+KXS/4)
             JX=J1-JXS*(KXS/4-KXT)
           CASE(3)
             IX=I1+IXS*(1+KXS/4-KXT)
             JX=J1+JXS*(1+KXS/4)
           CASE DEFAULT
             IX=I1-IXS*KXS/4
             JX=J1+JXS*(KXS/4-KXT)
           END SELECT
           NX=IJKGDS1(IX,JX,IJKGDSA)
           IF(NX.GT.0.)THEN
             IF(LI(NX,K).OR.IBI(K).EQ.0) THEN
               GO(N,K)=GI(NX,K)
               LO(N,K)=.TRUE.
               CYCLE N_LOOP
             ENDIF
           ENDIF
         ENDDO SPIRAL_LOOP
         IBO(K)=1
         GO(N,K)=0.
       ELSE
         IBO(K)=1
         GO(N,K)=0.
       ENDIF
     ELSE  ! no spiral search option
       IBO(K)=1
       GO(N,K)=0.
     ENDIF
   ENDDO N_LOOP
!$OMP END PARALLEL DO
 ENDDO KM_LOOP
 IF(IGDTNUMO.EQ.0) CALL POLFIXS(NO,MO,KM,RLAT,IBO,LO,GO)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 END SUBROUTINE POLATES3

max_stars_repo_name: noaa-ocs-modeling/PaHM

max_issues_repo_name: noaa-ocs-modeling/PaHM


text:
!> A container for the arrays that are used to store the distribution function among other things.
!!  These need to be accessible at a lower dependency level than the dist_fn module itself.
!! These arrays are allocated in the function dist_fn::allocate_arrays.

module dist_fn_arrays

   public :: gnew, gold, g_symm
   public :: g0, g1, g2, g3
   public :: g_krook, g_proj
   public :: gvmu
   public :: kperp2, dkperp2dr
   public :: wstar, wstarp
   public :: wdriftx_g, wdrifty_g
   public :: wdriftx_phi, wdrifty_phi
   public :: wdriftpx_g, wdriftpy_g
   public :: wdriftpx_phi, wdriftpy_phi

   ! dist fn
   complex, dimension(:, :, :, :, :), allocatable :: gnew, gold
   ! (naky, nakx, -nzgrid:nzgrid, ntubes, -vmu-layout-)

   complex, dimension(:, :, :, :, :), target, allocatable :: g_symm
   ! (naky, nakx, -nzgrid:nzgrid, ntubes, -vmu-layout-)

   complex, dimension(:, :, :, :, :), allocatable :: g0, g1, g2, g3
   ! (naky, nakx, -nzgrid:nzgrid, ntubes, -vmu-layout-)

   ! needed to implement time-delayed source when using Krook operator
   complex, dimension(:, :, :, :), allocatable :: g_krook
   ! (nakx, -nzgrid:nzgrid, ntubes, -vmu-layout-)

   ! needed to implement time-delayed source when using projection method
   complex, dimension(:, :, :, :), allocatable :: g_proj
   ! (nakx, -nzgrid:nzgrid, ntubes, -vmu-layout-)

   complex, dimension(:, :, :), allocatable :: gvmu
   ! (nvpa, nmu, -kxkyz-layout-)

   real, dimension(:, :, :), allocatable :: wstar, wstarp
   ! (nalpha, -nzgrid:nzgrid, -vmu-layout-)

   real, dimension(:, :, :), allocatable :: wdriftx_g, wdrifty_g
   real, dimension(:, :, :), allocatable :: wdriftx_phi, wdrifty_phi

   real, dimension(:, :, :), allocatable :: wdriftpx_g, wdriftpy_g
   real, dimension(:, :, :), allocatable :: wdriftpx_phi, wdriftpy_phi
   ! (nalpha, -nzgrid:nzgrid, -vmu-layout-)

   !> dkperp2dr will contain the radial variation of kperp2
   real, dimension(:, :, :, :), allocatable :: kperp2, dkperp2dr
   ! (naky, nakx, nalpha, -nzgrid:nzgrid)
   ! note: dkperp2dr is divided by kperp2

end module dist_fn_arrays

max_issues_repo_name: stellaGK/stella


text:
BLOCK DATA
END BLOCK DATA

max_stars_repo_name: ouankou/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
*> \brief \b SQPT01
*
*  =========== DOCUMENTATION ===========
*
* Online html documentation available at
*            http://www.netlib.org/lapack/explore-html/
*
*  Definition:
*  ===========
*
*       REAL             FUNCTION SQPT01( M, N, K, A, AF, LDA, TAU, JPVT,
*                        WORK, LWORK )
*
*       .. Scalar Arguments ..
*       INTEGER            K, LDA, LWORK, M, N
*       ..
*       .. Array Arguments ..
*       INTEGER            JPVT( * )
*       REAL               A( LDA, * ), AF( LDA, * ), TAU( * ),
*      $                   WORK( LWORK )
*       ..
*
*
*> \par Purpose:
*  =============
*>
*> \verbatim
*>
*> SQPT01 tests the QR-factorization with pivoting of a matrix A.  The
*> array AF contains the (possibly partial) QR-factorization of A, where
*> the upper triangle of AF(1:k,1:k) is a partial triangular factor,
*> the entries below the diagonal in the first k columns are the
*> Householder vectors, and the rest of AF contains a partially updated
*> matrix.
*>
*> This function returns ||A*P - Q*R||/(||norm(A)||*eps*M)
*> \endverbatim
*
*  Arguments:
*  ==========
*
*> \param[in] M
*> \verbatim
*>          M is INTEGER
*>          The number of rows of the matrices A and AF.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The number of columns of the matrices A and AF.
*> \endverbatim
*>
*> \param[in] K
*> \verbatim
*>          K is INTEGER
*>          The number of columns of AF that have been reduced
*>          to upper triangular form.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*>          A is REAL array, dimension (LDA, N)
*>          The original matrix A.
*> \endverbatim
*>
*> \param[in] AF
*> \verbatim
*>          AF is REAL array, dimension (LDA,N)
*>          The (possibly partial) output of SGEQPF.  The upper triangle
*>          of AF(1:k,1:k) is a partial triangular factor, the entries
*>          below the diagonal in the first k columns are the Householder
*>          vectors, and the rest of AF contains a partially updated
*>          matrix.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the arrays A and AF.
*> \endverbatim
*>
*> \param[in] TAU
*> \verbatim
*>          TAU is REAL array, dimension (K)
*>          Details of the Householder transformations as returned by
*>          SGEQPF.
*> \endverbatim
*>
*> \param[in] JPVT
*> \verbatim
*>          JPVT is INTEGER array, dimension (N)
*>          Pivot information as returned by SGEQPF.
*> \endverbatim
*>
*> \param[out] WORK
*> \verbatim
*>          WORK is REAL array, dimension (LWORK)
*> \endverbatim
*>
*> \param[in] LWORK
*> \verbatim
*>          LWORK is INTEGER
*>          The length of the array WORK.  LWORK >= M*N+N.
*> \endverbatim
*
*  Authors:
*  ========
*
*> \author Univ. of Tennessee
*> \author Univ. of California Berkeley
*> \author Univ. of Colorado Denver
*> \author NAG Ltd.
*
*> \date November 2011
*
*> \ingroup single_lin
*
*  =====================================================================
      REAL             FUNCTION SQPT01( M, N, K, A, AF, LDA, TAU, JPVT,
     $                 WORK, LWORK )
      implicit none
*
*  -- LAPACK test routine (version 3.4.0) --
*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
*     November 2011
*
*     .. Scalar Arguments ..
      INTEGER            K, LDA, LWORK, M, N
*     ..
*     .. Array Arguments ..
      INTEGER            JPVT( * )
      REAL               A( LDA, * ), AF( LDA, * ), TAU( * ),
     $                   WORK( LWORK )
*     ..
*
*  =====================================================================
*
*     .. Parameters ..
      REAL               ZERO, ONE
      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
*     ..
*     .. Local Scalars ..
      INTEGER            I, INFO, J
      REAL               NORMA
*     ..
*     .. Local Arrays ..
      REAL               RWORK( 1 )
*     ..
*     .. External Functions ..
      REAL               SLAMCH, SLANGE
      EXTERNAL           SLAMCH, SLANGE
*     ..
*     .. External Subroutines ..
      EXTERNAL           SAXPY, SCOPY, SORMQR, XERBLA
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          MAX, MIN, REAL
*     ..
*     .. Executable Statements ..
*
      SQPT01 = ZERO
*
*     Test if there is enough workspace
*
      IF( LWORK.LT.M*N+N ) THEN
         CALL XERBLA( 'SQPT01', 10 )
         RETURN
      END IF
*
*     Quick return if possible
*
      IF( M.LE.0 .OR. N.LE.0 )
     $   RETURN
*
      NORMA = SLANGE( 'One-norm', M, N, A, LDA, RWORK )
*
      DO 30 J = 1, K
         DO 10 I = 1, MIN( J, M )
            WORK( ( J-1 )*M+I ) = AF( I, J )
   10    CONTINUE
         DO 20 I = J + 1, M
            WORK( ( J-1 )*M+I ) = ZERO
   20    CONTINUE
   30 CONTINUE
      DO 40 J = K + 1, N
         CALL SCOPY( M, AF( 1, J ), 1, WORK( ( J-1 )*M+1 ), 1 )
   40 CONTINUE
*
      CALL SORMQR( 'Left', 'No transpose', M, N, K, AF, LDA, TAU, WORK,
     $             M, WORK( M*N+1 ), LWORK-M*N, INFO )
*
      DO 50 J = 1, N
*
*        Compare i-th column of QR and jpvt(i)-th column of A
*
         CALL SAXPY( M, -ONE, A( 1, JPVT( J ) ), 1, WORK( ( J-1 )*M+1 ),
     $               1 )
   50 CONTINUE
*
      SQPT01 = SLANGE( 'One-norm', M, N, WORK, M, RWORK ) /
     $         ( REAL( MAX( M, N ) )*SLAMCH( 'Epsilon' ) )
      IF( NORMA.NE.ZERO )
     $   SQPT01 = SQPT01 / NORMA
*
      RETURN
*
*     End of SQPT01
*
      END


text:
      subroutine wshrun
c
c     + + + PURPOSE + + +
c
c     SR WSHRUN routes the runoff from impoundments and
c     channels through the watershed.
c
c     Called from: SR WSHDRV
c     Author(s): Ascough II
c     Reference in User Guide:
c
c     Version:
c     Date recoded:
c     Recoded by:
c
c     + + + KEYWORDS + + +
c
c     + + + PARAMETERS + + +
c
      include 'pmxcsg.inc'
      include 'pmxelm.inc'
      include 'pmxhil.inc'
      include 'pmximp.inc'
      include 'pmxpln.inc'
      include 'pmxprt.inc'
c
c     + + + ARGUMENT DECLARATIONS + + +
c
c     + + + ARGUMENT DEFINITIONS + + +
c
c     + + + COMMON BLOCKS + + +
c
      include 'cdiss11.inc'
      include 'cirspri.inc'
      include 'coutchn.inc'
      include 'cstore1.inc'
      include 'cstruct.inc'

      include 'chydrol.inc'
      include 'cchpek.inc'
      include 'cstruc.inc'
c
c     + + + LOCAL VARIABLES + + +
c
c     + + + LOCAL DEFINITIONS + + +
c
c     + + + SAVES + + +
c
c     + + + SUBROUTINES CALLED + + +
c
c     wshpek
c
c     + + + DATA INITIALIZATIONS + + +
c
c     + + + END SPECIFICATIONS + + +
c
c
c     calculate the average runoff on the watershed at the outlet
c     of the current channel element
c
      rofave(ielmt) = runvol(ielmt) / wsarea(ielmt)
c
c     save monthly, annual, and total runoff from each channel
c     element
c
      trunm(ielmt) = trunm(ielmt) + runvol(ielmt)
      truny(ielmt) = truny(ielmt) + runvol(ielmt)
      trunt(ielmt) = trunt(ielmt) + runvol(ielmt)
cd    Added by S. Dun, Dec 20, 2007
c     for more output of annual summary of the chnnel
      tronvy(ielmt) = tronvy(ielmt) + rvolon(ielmt)
      tronvt(ielmt) = tronvt(ielmt) + rvolon(ielmt)
      sbrvty(ielmt) = sbrvty(ielmt) + sbrunf(ielmt-nhill)
     1               * charea(ielmt-nhill)
      sbrvtt(ielmt) = sbrvtt(ielmt) + sbrunf(ielmt-nhill)
     1               * charea(ielmt-nhill)
c     End adding
c
c     the duration for the channel element is the maximum
c     duration of the watershed elements that may feed it
c
cd    Modified by S. Dun 1/08/2004
cd    This duration value has given in WSHCQI already.
cd    Only the dur for channel could be changed after WSHIRS.for
cd
cd      watdur(ielmt) = amax1(watdur(nhleft(ielmt)),
cd     1    watdur(nhrght(ielmt)),watdur(nhtop(ielmt)),
cd     1    watdur(ncleft(ielmt)),watdur(ncrght(ielmt)),
cd     1    watdur(nctop(ielmt)),watdur(nileft(ielmt)),
cd     1    watdur(nirght(ielmt)),watdur(nitop(ielmt)),dur,irdur)
      watdur(ielmt) = amax1(watdur(ielmt),dur)
cd    End modifying
c
c
c     calculate peak runoff using modified Rational method
c     (ipeak = 1) or modified CREAMS method (ipeak = 2)
c
      call wshpek
c
      return
      end

ext: for

max_stars_repo_name: jarad/dep

max_issues_repo_name: jarad/dep

max_forks_repo_name: akrherz/idep


text:
MODULE CRTM_SharedData

! Provides CRTM-wide access to shared data without the
! need to USE specific shared data modules.

  USE CRTM_SpcCoeff    , ONLY: SC
  USE CRTM_TauCoeff    , ONLY: TC
  USE CRTM_CloudCoeff  , ONLY: CloudC
  USE CRTM_AerosolCoeff, ONLY: AeroC
  USE CRTM_EmisCoeff   , ONLY: EmisC
  USE CRTM_BeCoeff     , ONLY: BeC

END MODULE CRTM_SharedData

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


text:
C     This program is written for testing module TYPING
C
C
      PROGRAM TYPAGE_INTRINSIC
      LOGICAL L
      INTEGER I
      REAL*4 R
      REAL*8 D
      COMPLEX*8 C
      COMPLEX*16 E
      CHARACTER CH

C     Conversion intrinsics
C     =====================
      I = INT(R)
      I = INT(D)
      I = INT(C)
      I = INT(CH)               ! ERROR
      I = INT(L)                ! ERROR

      R = REAL(I)
      R = REAL(D)
      R = REAL(C)
      R = REAL(CH)              ! ERROR
      R = REAL(L)               ! ERROR

      D = DBLE(I)
      D = DBLE(R)
      D = DBLE(C)
      D = DBLE(CH)              ! ERROR
      D = DBLE(L)               ! ERROR

      C = CMPLX(I)
      C = CMPLX(R)
      C = CMPLX(D)
      C = CMPLX(CH)             ! ERROR
      C = CMPLX(L)              ! ERROR

C     Argument & Result: REAL or DOUBLE
C     ========================
      I = TAN(I)                ! ERROR
      R = TAN(R)
      D = TAN(D)
      I = TAN(C)                ! ERROR
      I = TAN(L)                ! ERROR

C     Argument & Result: REAL, DOUBLE or COMPLEX
C     =================================
      I = EXP(I)                ! ERROR
      R = EXP(R)
      D = EXP(D)
      C = EXP(C)
      I = EXP(L)                ! ERROR

      I = LOG(I)                ! ERROR
      R = LOG(R)
      D = LOG(D)
      C = LOG(C)
      I = LOG(L)                ! ERROR

C     Argument & Result: INT, REAL or DOUBLE
C     =============================
      I = MIN(I, 15, I)
      R = MIN(I, 15, R)
      D = MIN(I, 15, D)
      C = MIN(I, 15, C)         ! ERROR
      I = MIN(I, 15, L)         ! ERROR

      I = MOD(I, I)
      R = MOD(R, R)
      D = MOD(D, D)
      I = MOD(R, D)
      C = MOD(C, C)             ! ERROR

C     Intrinsic ABS: IRDC -> IRDR
C     ===========================
      I = ABS(I)
      R = ABS(R)
      D = ABS(D)
      R = ABS(C)
      C = ABS(C)
      I = ABS(L)                ! ERROR

C     Type_argument --> Type_result
C     ===================================
C     INT  ->  INT
      I = IABS(I)
      I = IABS(R)               ! ERROR
      I = IABS(D)               ! ERROR
      I = IABS(C)               ! ERROR

C     REAL -> REAL
      R = ALOG(I)               ! ERROR
      R = ALOG(R)
      R = ALOG(D)               ! ERROR
      R = ALOG(C)               ! ERROR

C     DOUBLE -> DOUBLE
      D = DLOG(I)               ! ERROR
      D = DLOG(R)               ! ERROR
      D = DLOG(D)
      D = DLOG(C)               ! ERROR

C     COMPLEX -> COMPLEX
      C = CLOG(I)               ! ERROR
      C = CLOG(R)               ! ERROR
      C = CLOG(D)               ! ERROR
      C = CLOG(C)

C     CHAR -> INT
      I = ICHAR(CH)
      D = ICHAR(R)              ! ERROR

      I = LEN(CH)
      I = LEN(R)                ! ERROR

      I = INDEX(CH)
      I = INDEX(R)              ! ERROR

C     INT      ->   CHAR
      CH = CHAR(I)
      CH = CHAR(D)              ! ERROR

C     REAL     ->   INT
      I = MIN1(I, I)            ! ERROR
      I = MIN1(R, R)
      I = MIN1(I, R)            ! ERROR
      I = MIN1(I, R, D)         ! ERROR
      I = MIN1(I, R, D, C)      ! ERROR

C     COMPLEX  ->   REAL
      R = AIMAG(I)              ! ERROR
      R = AIMAG(R)              ! ERROR
      R = AIMAG(D)              ! ERROR
      R = AIMAG(C)

C     REAL     ->   DOUBLE
      D = DPROD(I)              ! ERROR
      D = DPROD(R)
      D = DPROD(D)              ! ERROR
      D = DPROD(C)              ! ERROR

C     INT      ->   REAL
      R = AMIN0(I, I)
      R = AMIN0(I, R)           ! ERROR
      R = AMIN0(R, D)           ! ERROR
      R = AMIN0(D, C)           ! ERROR

C     CHAR     ->   LOGICAL
      L = LGE(I, CH)            ! ERROR
      L = LGE(R, CH)            ! ERROR
      L = LGE(CH, CH)

      END


text:
C Copyright (C) 2002, Carnegie Mellon University and others.
C All Rights Reserved.
C This code is published under the Common Public License.

      subroutine FACT_FW(TASK, NZ, NY, NCOL, NAC, LA48, NNZ, A, IRA,
     1     JCA, A48, IRN48, JCN48, KEEP48, LRW, RW, LIW, IW, IERR)
C
C     $Id: fact_fw.f,v 1.1 2002/05/02 18:52:15 andreasw Exp $
C
C     Factorize basis matrix for one finite element
C
C     TASK   1: call MA48AD to prepare datastructures
C               in:  LA48, NNZ, A(NNZ), IRA(NNZ), JCA(NNZ) (all unchanged)
C               out: IRN48(LA48), JCN48(NNZ), KEEP48(*)
C            2: call MA48BD with job = 1 (i.e. KEEP and IRN can change)
C               in:  LA48, NNZ, A(NNZ), IRN48(LA48), JCN48(NNZ), KEEP48()
C               out: A48(LA48), IRN48(LA48), KEEP48()
C            3: call MA48BD with job = 2 (fast call), doesn't change KEEP, IRN
C               in:  LA48, NNZ, A(NNZ), IRN48(LA48), JCN48(NNZ), KEEP48()
C               out: A48(LA48)
C
C     Author:  Andreas Waechter
C              c/o Group of Larry Biegler
C              Department of Chemical Engineering
C              Carnegie Mellon University
C              Pittsburgh, PA
C
      implicit none

      integer TASK, NZ, NY, NCOL, NAC, LA48, NNZ
      double precision A(NNZ)
      integer IRA(NNZ)
      integer JCA(NNZ)
      double precision A48(LA48)
      integer IRN48(LA48)
      integer JCN48(NNZ)
      integer KEEP48(7+10*NCOL*(NZ+NY))
      integer LRW
      double precision RW(LRW)
      integer LIW
      integer IW(LIW)
      integer IERR

      integer p_rwend, p_iwend, p_a, p_iw, p_w, p_jcn
      integer ntot, job, i
      integer info(12), icntl(9)
      double precision rinfo, cntl(5)

      p_rwend = 0
      p_iwend = 0
      IERR    = 0
      call MA48ID(cntl, icntl)
      icntl(2) = 0
C
C     The following statement changes the pivot tolerance
C
c      cntl(2)  = 1.d-2

      goto (1000, 2000, 3000) TASK

c      write(*,*) 'wrong TASK = ',TASK
      STOP
C
C     prepare data structures for MA48AD
C
 1000 continue

      job = 1
      ntot = NCOL*(NZ+NY) + NAC

      p_a     = p_rwend
      p_rwend = p_a + LA48
      if( p_rwend.gt.LRW ) goto 9098
      call DCOPY(NNZ, A, 1, RW(p_a+1), 1)

      p_iw    = p_iwend
      p_jcn   = p_iw  + 9*ntot
      p_iwend = p_jcn + LA48
      if( p_iwend.gt.LIW ) goto 9099

      do i = 1, NNZ
         IRN48(i)    = IRA(i)
         IW(p_jcn+i) = JCA(i)
      enddo

      call MA48AD(ntot, ntot, NNZ, job, LA48, RW(p_a+1), IRN48,
     1     IW(p_jcn+1), KEEP48, cntl, icntl, IW(p_iw+1), info, rinfo)
      if( info(1).ne.0 ) then
         if( info(1).eq.-4 .or. info(1).eq.2 ) then
            IERR = 20           ! system singular
            goto 9999
         endif
         if( info(1).eq.-3 ) then
            IERR = 21           ! LA48 too small
            goto 9999
         endif
         if( info(1).ne.1 .or. info(12).ne.0 ) then
c            write(*,*) 'info = ',info(1)
            IERR = 48
            goto 9999
         endif
      endif

      do i = 1, NNZ
         JCN48(i) = IW(p_jcn+i)
      enddo

      p_rwend = p_a
      p_iwend = p_iw

      goto 9999
C
C     Do normal factorization
C
 2000 continue
      job = 1
      goto 3500
C
C     Do fast factorization
C
 3000 continue
      job = 2

 3500 continue
      ntot = NCOL*(NZ+NY) + NAC

      p_w     = p_rwend
      p_rwend = p_w + ntot
      if( p_rwend.gt.LRW ) goto 9098

      p_iw    = p_iwend
      p_iwend = p_iw + 4*ntot
      if( p_iwend.gt.LIW ) goto 9099

      call DCOPY(NNZ, A, 1, A48, 1)

      call MA48BD(ntot, ntot, NNZ, job, LA48, A48, IRN48, JCN48, KEEP48,
     1     cntl, icntl, RW(p_w+1), IW(p_iw+1), info, rinfo)
      if( info(1).ne.0 ) then
         if( info(1).eq.-3 ) then
            IERR = 21
            goto 9999
         endif
c         write(*,*) 'info = ',info(1)
         IERR = 48
         goto 9999
      endif

      goto 9999
C
C     Error handling
C
 9098 IERR = 98
      goto 9999
 9099 IERR = 99
      goto 9999

 9999 continue
      return
      end

max_stars_repo_name: kant/CMMPPT

max_issues_repo_name: kant/CMMPPT

max_forks_repo_name: kant/CMMPPT


text:
As of January 2007, Yolo County Superior Court will accept online traffic school.  Make sure to check with the county you received your citation from. You can select an online traffic school from their list.

I took traffic school online.  It was a few years ago, but I think I just Googled Traffic School Yolo County or something like that.  It took a couple hours. Users/JanelleAlvstadMattson

Ive done both traffic schools, online and in person classes.  Online took me a couple hours too and in person took one Saturday from 8 to 5. Both cost around the same amount around 30 to 40 bucks, just one takes a Saturday.  But I liked online more because I got to do it on my time. For Solano County if you do traffic school online you have to take a test conducted by the NTSA in addition to your online traffic school of choice.  The final test is really easy if you can flip through the California drivers manual and find anwers.  Its exactly the type of questions from a written exam to get a license. Lori Tobias

Ive used http://www.ampmtrafficschool.com AMPMTrafficSchool and I was really happy with it.  Quick, easy, and cheap.  Users/PatrickProctor

20060515 17:10:40 nbsp I took Woodland Adult Education.  8:30 to 4:00 with a half hour lunch break.  The teacher was a big, surly, excop, but he was pretty interesting.  Cost was $30.  Make sure you sign up at least a week in advance.  They get a lot more business now that you cant go online. Users/VihilaAlBhed

20070117 14:17:52 nbsp I received all the forms from the Yolo County Superior court to attend a traffic school. As of December 2006, you can attend an online traffic school, again. The line saying The Yolo SC... does NOT allow home study traffic school... is striked through. Users/GianlucaVarenni

20070129 20:08:05 nbsp Go to Cheap School, Ada Vides rocks! Shes a hilarious New York sounding lady, kept us pretty entertained and sO Nice. You dont even have to take a real quiz! Users/TheStomach

20080229 15:33:27 nbsp Just make sure that the traffic school that you choose is a court approved one or you wont receive credit. Make sure that the court gives you a list (theyre supposed to). I did AMPM a while back and it was cake.   Users/CurlyGirl26

20080313 16:08:28 nbsp The information on this page is wrong. The Yolo County reversed its policy and it accepts online traffic school. I recommend ampmtrafficschool.com They are local and provide great service. Users/Safa

20090211 11:52:25 nbsp I used ampm, which was on the list approved by the court. Was easy, quick and cheap ($17.99) Users/DanielZinn

20100129 12:43:16 nbsp Question...so in Yolo you have to pay the $66 to essentailly be granted the ability to take traffic school, and then you also have to pay for the online course which is seperate from the $66 fee, correct?  Users/AmLin
  Thats correct, but generally its still much cheaper than the change in your insurance premium would be.  Users/rfrazier
  Will traffic school essentially erase the ticket? I know it takes the point back off your record, but if you do the traffic school will your insurance remain unaffected? Any experience with this? AmLin
  I have taken traffic school a couple of times and yes, it does keep it off your record.  Essentially the ticket is erased, however keep in mind that you can only do it every 18 months.  RF
   You can sometimes also take extended traffic school for a singlepoint violation even within that 18 month period.  It costs still more, though.  Traffic school will clear your record of the citation.  Also note that a single point on your record shouldnt increase your rates.  Onepoint and some twopoint violations are cleared from your record without traffic school after 3 years.  Some offenses such as DUI, or lesser offenses exacerbated by failure to show up at court will last longer5 or 10 years, depending. Users/TomGarberson
   It was my first ticket ever, so I will definitely be trying to stay out of trouble for the next 18 months! Thanks for the info guys. AmLin

20100412 12:31:01 nbsp Based on the comments below I took AMPMtraffic online school. It was on the list of schools from the Yolo Court. Took about 3 hours. They include offsubject questions within the material like Robert Anderson of Scotland invented the first crude electric car, which was not named after him to make sure you arent gaming the school and not reading the material. Remember there there are ten classes of licenses and one question involves 42% versus 15 % on accidents by drunk/sober drivers. You have to get 100% on the quizzes and 80% on the final.  Cost was just under $20, and I dont yet know if the court got/processed the certificate. My dealings with the court have been nothing less than frustrating. Due to the cuts in government budgets, the court has one furlough day a month, the ticket wasnt registered so I could apply for traffic school until I phoned twice and emailed once. My advice is to stay out of the eye of the law as I almost had to go to court even though I wanted traffic school due to these problems.  Users/ArgusDemosthenes


text:
module drivers

! This module contains high level drivers for atomic SCF calculations
!
! The idea is to use these drivers to do most frequent calculations with
! an exponential mesh and to get an idea how things work. They can be used as a
! starting point/template to write a custom solver for a particular problem,
! or to use a different mesh.

use types, only: dp
use constants, only: pi
use ode1d, only: integrate
use states, only: get_atomic_states_nonrel, get_atomic_states_rel
use dft, only: KS_step, rho2V
use mesh, only: mesh_exp, mesh_exp_deriv
use energies, only: get_hydrogen_energies, thomas_fermi_potential, E_nl
use dft_data, only: dft_data_t
use mixings, only: mixing_anderson
implicit none

private
public get_atom_orb, atom_lda, atom_rlda, get_atom_orb_rel

contains

integer function get_atom_orb(Z) result(n_orb)
! Returns the number of orbitals for the given atom
integer, intent(in) :: Z

integer, pointer, dimension(:) :: no, lo
real(dp), pointer, dimension(:) :: fo

call get_atomic_states_nonrel(Z, no, lo, fo)
n_orb = size(no)
deallocate(no, lo, fo)
end function

integer function get_atom_orb_rel(Z) result(n_orb)
! Returns the number of orbitals for the given atom
integer, intent(in) :: Z

integer, pointer, dimension(:) :: no, lo, so
real(dp), pointer, dimension(:) :: fo

call get_atomic_states_rel(Z, no, lo, so, fo)
n_orb = size(no)
deallocate(no, lo, so, fo)
end function

subroutine atom_lda(Z, r_min, r_max, a, N, no, lo, fo, ks_energies, E_tot, &
    R, Rp, V_tot, density, orbitals, reigen_eps, reigen_max_iter, &
    mixing_eps, mixing_alpha, &
    mixing_max_iter, perturb)
integer, intent(in) :: Z
real(dp), intent(in) :: r_min, r_max, a
integer, intent(in) :: N
integer, intent(out), target :: no(:), lo(:)
real(dp), intent(out), target :: fo(:)
real(dp), intent(out), target :: ks_energies(:)
real(dp), intent(out) :: E_tot
real(dp), intent(out), target :: R(:), Rp(:)
real(dp), intent(out), target :: V_tot(:)
real(dp), intent(out), target :: density(:)
real(dp), intent(out), target :: orbitals(:, :)
real(dp), intent(in) :: reigen_eps, mixing_eps, mixing_alpha
integer, intent(in) :: mixing_max_iter, reigen_max_iter
logical, intent(in) :: perturb

integer, pointer :: no_a(:), lo_a(:)
real(dp), pointer :: fo_a(:)
integer :: i
real(dp), dimension(size(ks_energies)), target :: Emin_init, Emax_init

real(dp), dimension(size(R)), target :: V_h, V_xc, e_xc, V_coulomb, tmp
type(dft_data_t) :: d

call get_atomic_states_nonrel(Z, no_a, lo_a, fo_a)
no = no_a
lo = lo_a
fo = fo_a
deallocate(no_a, lo_a, fo_a)

R = mesh_exp(r_min, r_max, a, N)
Rp = mesh_exp_deriv(r_min, r_max, a, N)

!V_tot = -Z / R
ks_energies = get_hydrogen_energies(Z, no)
V_tot = thomas_fermi_potential(R, Z)
!ks_energies = get_tf_energies(Z, no, fo)

V_coulomb = -Z/R

! We allow a few unbounded states
Emax_init = 10
! Use Hydrogen Schroedinger energies for the lower limit:
do i = 1, size(no)
    Emin_init(i) = E_nl(0._dp, no(i), lo(i), Z, 0)
end do
! For robustness, decrease Emin by 10%:
Emin_init = 1.1_dp * Emin_init

d%Z = Z
d%R => R
d%Rp => Rp
d%rho => density
d%V_h => V_h
d%V_coulomb => V_coulomb
d%V_xc => V_xc
d%e_xc => e_xc
d%V_tot => V_tot
d%orbitals => orbitals
d%alpha = mixing_alpha
d%dirac = .false.
d%perturb = perturb
d%reigen_eps = reigen_eps
d%reigen_max_iter = reigen_max_iter
d%no => no
d%lo => lo
d%fo => fo
d%ks_energies => ks_energies
d%Emax_init => Emax_init
d%Emin_init => Emin_init
! Start from an initial density instead:
!d%rho = 1 / cosh(d%R)**2
!d%rho = d%rho * d%Z / integrate(d%Rp, 4*pi*d%rho*d%R**2)
!call rho2V(d)
tmp = mixing_anderson(KS_step, d%V_tot - d%V_coulomb, mixing_max_iter, &
    .true., d, mixing_alpha, mixing_eps)
E_tot = d%Etot
! Prints the energies:
!print *, "Ekin: ", d%Ekin
!print *, "Ecoul:", d%Ecoul
!print *, "Eenuc:", d%Eenuc
!print *, "Exc:  ", d%Exc
!print *, "E_tot:", d%Ekin + d%Ecoul + d%Eenuc + d%Exc
end subroutine

subroutine atom_rlda(Z, r_min, r_max, a, N, c, no, lo, so, fo, ks_energies, &
    E_tot, &
    R, Rp, V_tot, density, orbitals, reigen_eps, reigen_max_iter, &
    mixing_eps, mixing_alpha, &
    mixing_max_iter, perturb)
integer, intent(in) :: Z
real(dp), intent(in) :: r_min, r_max, a, c
integer, intent(in) :: N
integer, intent(out), target :: no(:), lo(:), so(:)
real(dp), intent(out), target :: fo(:)
real(dp), intent(out), target :: ks_energies(:)
real(dp), intent(out) :: E_tot
real(dp), intent(out), target :: R(:), Rp(:)
real(dp), intent(out), target :: V_tot(:)
real(dp), intent(out), target :: density(:)
real(dp), intent(out), target :: orbitals(:, :)
real(dp), intent(in) :: reigen_eps, mixing_eps, mixing_alpha
integer, intent(in) :: mixing_max_iter, reigen_max_iter
logical, intent(in) :: perturb

integer, pointer :: no_a(:), lo_a(:), so_a(:)
real(dp), pointer :: fo_a(:)
real(dp), dimension(size(ks_energies)), target :: Emin_init, Emax_init
integer :: i, relat

real(dp), dimension(size(R)), target :: V_h, V_xc, e_xc, V_coulomb, tmp
type(dft_data_t) :: d

call get_atomic_states_rel(Z, no_a, lo_a, so_a, fo_a)
no = no_a
lo = lo_a
so = so_a
fo = fo_a
deallocate(no_a, lo_a, so_a, fo_a)

R = mesh_exp(r_min, r_max, a, N)
Rp = mesh_exp_deriv(r_min, r_max, a, N)

!V_tot = -Z / R
ks_energies = get_hydrogen_energies(Z, no)
V_tot = thomas_fermi_potential(R, Z)
!ks_energies = get_tf_energies(Z, no, fo)

V_coulomb = -Z/R

! We allow a few unbounded states
Emax_init = 10
! Use Hydrogen Dirac energies for the lower limit:
do i = 1, size(no)
    if (so(i) == 1) then
        relat = 2
    else
        relat = 3
    end if
    Emin_init(i) = E_nl(c, no(i), lo(i), Z, relat)
end do
! For robustness, decrease Emin by 10%:
Emin_init = 1.1_dp * Emin_init

d%Z = Z
d%R => R
d%Rp => Rp
d%rho => density
d%V_h => V_h
d%V_coulomb => V_coulomb
d%V_xc => V_xc
d%e_xc => e_xc
d%V_tot => V_tot
d%orbitals => orbitals
d%alpha = mixing_alpha
d%dirac = .true.
d%perturb = perturb
d%reigen_eps = reigen_eps
d%reigen_max_iter = reigen_max_iter
d%no => no
d%lo => lo
d%fo => fo
d%ks_energies => ks_energies
d%Emax_init => Emax_init
d%Emin_init => Emin_init

! Relativistic parameters:
d%c = c
d%so => so

tmp = mixing_anderson(KS_step, d%V_tot - d%V_coulomb, mixing_max_iter, &
    .true., d, mixing_alpha, mixing_eps)
E_tot = d%Etot
end subroutine

end module

max_stars_repo_name: rc/dftatom

max_issues_repo_name: rc/dftatom

max_forks_repo_name: rc/dftatom


text:
module proj_parameters

   implicit none

   integer, parameter :: initial_projection_comp = 1
   integer, parameter :: divu_iters_comp         = 2
   integer, parameter :: pressure_iters_comp     = 3
   integer, parameter :: regular_timestep_comp   = 4

end module proj_parameters

max_stars_repo_name: sailoridy/MAESTRO

max_issues_repo_name: sailoridy/MAESTRO

max_forks_repo_name: sailoridy/MAESTRO