text:
c=======================================================================
c
c     subroutine RETEM
c
c     Expectation-Maximization (EM) algorithm
c     -> recover missing observations in a time series
c
c     Return data structure
c
c           n : number of points
c           p : number of assets
c
c       R(n,p) = | x(1,1), ..., x(1,p) |
c                | x(2,1), ..., -1000  |
c                | x(3,1), ..., x(3,p) |
c                |   .   , ...,   .    |
c                |   .   , ...,   .    |
c                | x(n,1), ..., x(n,p) |
c
c     where x(2,p) <= -1000 denote a missing return
c
c-----------------------------------------------------------------------
      SUBROUTINE retem ( n, p, X, miss, iwork, dwork, Y, info )
c-----------------------------------------------------------------------
c
c     INPUT
c        n      : max. number of values (n > 1)                  integer
c        p      : number of assets (p > 0)                       integer
c        X      : returns values (n*p)                            double
c        miss   : missing value (<= -1)                           double
c
c     WORKSPACE
c        iwork  : 2*p                                            integer
c        dwork  : p*(13*p + 2*n + 10)                             double
c
c     OUTPUT
c        Y      : EM-recoved returns values (n*p)                 double
c        info   : diagnostic argument                            integer
c
c     CALL
c
c-----------------------------------------------------------------------
c
      IMPLICIT NONE
c
c     arguments i/o
      INTEGER n, p, info
      DOUBLE PRECISION miss, X(n,*), Y(n,*)
c
c     workspaces
      INTEGER iwork(*)
      DOUBLE PRECISION dwork(*)
c
c     local variables
      INTEGER t, i, j, k, l, m, nobs, mt, nt, mtt, h
      INTEGER MaxLoop, CountLoop, Convergence
      INTEGER pindMiss, pindObs
      INTEGER pdX, pdy,
     &        pdrho, pdrhoN, pdrhoMiss, pdrhoObs,
     &        pdcov, pdcovN, pdcovm,
     &        pdyMiss, pdyObs, pdM, pdyM, pdv,
     &        pdcovMO, pdcovMM, pdcovOO, pdinvCovOO, pdjms,
     &        pdCC, pdM1, pdM2, pdM3, pdMt, pdC
      DOUBLE PRECISION Tolerance, Distance, tmp, EPS, ZERO
      PARAMETER (EPS = 1.E-15, ZERO = 0.0)
c
c     external subroutines
      EXTERNAL IVX, YM, MCM, covm
c
c-----------------------------------------------------------------------
c
c     initialization
      info = 0
      CALL YM ( n, p, X, Y ) ! copy X -> Y
c
c     pointers for integer work space  : iwork
c     ----------------------------------------
      pindMiss = 1
c     pindMiss : pointer missing data index (p)
      pindObs  = pindMiss + ( p )
c     pindObs  : pointer observed data index (p)
c
c     Total size of iwork array  = 2*p
c
c     pointers for double precision work space  : dwork
c     -------------------------------------------------
      pdX    = 1
c     pdX    : pointer for availaible NAV (n*p)
      pdrho  = pdX + ( n*p )
c     pdrho  : pointer for mean vector (p)
      pdrhoN = pdrho + ( p )
c     pdrhoN : pointer for new mean vector (p)
      pdcov  = pdrhoN + ( p )
c     pdcov  : pointer for covariance matrix (p*p)
      pdcovN = pdcov + ( p*p )
c     pdcovN : pointer for new covariance matrix (p*p)
      pdcovm = pdcovN + ( p*p )
c     pdcovm : pointer for COVM (p*n + p)
      pdy    = pdcovm + ( n*p + p)
c     pdy    : pointer (p)
      pdyMiss = pdy + ( p )
c     pdyMiss : pointer missing data y (p)
      pdyObs  = pdyMiss + ( p )
c     pdyObs  : pointer observed data y (p)
      pdM     = pdyObs + ( p )
c
      pdyM    = pdM + ( p )
c
      pdv     = pdyM + ( p )
c
      pdrhoMiss  = pdv + ( p )
c     pdrhoMiss  : pointer missing mean (p)
      pdrhoObs   = pdrhoMiss + ( p )
c     pdrhoObs   : pointer observed mean (p)
      pdcovMO    = pdrhoObs + ( p )
c     pdcovMO    : pointer missing/observed cov (p*p)
      pdcovMM    = pdcovMO + ( p*p )
c     pdcovMM    : pointer missing/missing cov (p*p)
      pdcovOO    = pdcovMM + ( p*p )
c     pdcovOO    : pointer observed/observed cov (p*p)
      pdinvCovOO = pdcovOO + ( p*p )
c     pdinvCovOO : pointer observed/observed inverse cov (p*p)
      pdCC       = pdinvCovOO + ( p*p )
c
      pdM1       = pdCC + ( p*p )
c
      pdM2       = pdM1 + ( p*p )
c
      pdM3       = pdM2 + ( p*p )
c
      pdMt       = pdM3 + ( p*p )
c
      pdC        = pdMt + ( p*p )
c
      pdjms      = pdC + ( p*p )
c     pdjms      : pointer for JMS (p*p)
c
c     Total size of dwork array  = p*(13*p + 2*n + 10)
c
c-----------------------------------------------------------------------
c
c     E-M initialization
      h  = 1       ! horizon (log-returns
      m  = 0       ! number of missing data (one by row)
      DO t = 1,n
        DO i = 1,p
            IF (X(t,i) .LE. miss) THEN
                m = m + 1
                GOTO 100
            ENDIF
        ENDDO
  100 CONTINUE
      ENDDO
      IF (m .EQ. 0) RETURN ! if no missing data -> exit
      nobs = n - m         ! number of available dat
      IF (nobs .EQ. 0) THEN
        info = -9003 ! no enough available data for initialization
        RETURN
      ENDIF
c
      k = 0
      DO t = 1,n
        DO i = 1,p
            IF (X(t,i) .LE. miss) GOTO 101
        ENDDO
        k = k + 1
        DO i = 1,p
            dwork(pdX + nobs*(i-1) + k - 1) = X(t,i)
        ENDDO
  101 CONTINUE
      ENDDO
c
c     mean and cov. matrix
      CALL MCM (nobs,p,dwork(pdX),dwork(pdrho))
      CALL YV (p,dwork(pdrho),dwork(pdrhoN))
      CALL covm (nobs,p,dwork(pdX),dwork(pdcovm),dwork(pdcov),info)
c
      CALL YV (p*p,dwork(pdcov),dwork(pdcovN))
      IF (info .LT. 0) RETURN
c
c     EM loop
      tmp = dwork(pdrhoN) + SQRT(MAX(dwork(pdcovN),EPS))
      DO i = 2,p
        tmp = tmp + dwork(pdrhoN+i-1) +
     &        SQRT(MAX(dwork(pdcovN+p*(i-1)+i-1),EPS))
      ENDDO
      Tolerance = 1.E-6*(tmp/(2.0*p)) ! algorithm tolerance
      MaxLoop     = 100  ! maximum EM loop
      CountLoop   = 0    ! loop counter
      Convergence = 0    ! convergence test (0 no, 1 yes)
      DO WHILE ((Convergence .EQ. 0).AND.(CountLoop.LT.MaxLoop))
        CountLoop = CountLoop + 1  ! counter + 1
        CALL YV (p, dwork(pdrhoN), dwork(pdrho))
        CALL YV (p*p, dwork(pdcovN), dwork(pdcov))
c
c       Step 1: estimation
        CALL IVX (p*p, dwork(pdC), ZERO)
        CALL IVX (p*p, dwork(pdMt), ZERO) ! Mt = zeros(p,p)
        DO t = 1,n
             mt = 0 ! number of missing data at time t
             nt = 0 ! number of observed data at time t
             CALL IVX (p*p, dwork(pdCC), ZERO) ! CC = zeros(p,p)
             DO i = 1,p
                dwork(pdy + i - 1) = X(t,i)
                IF (X(t,i) .LE. miss) THEN
                    dwork(pdyMiss + mt)   = dwork(pdy+i-1)
                    iwork(pindMiss + mt)  = i ! index of missing data
                    dwork(pdrhoMiss + mt) = dwork(pdrho + i - 1)
                    mt = mt + 1 ! number of missing data at time t
                ELSE
                    dwork(pdyObs + nt)   = dwork(pdy+i-1) ! observation
                    iwork(pindObs + nt)  = i ! index of observed data
                    dwork(pdrhoObs + nt) = dwork(pdrho + i - 1)
                    nt = nt + 1 ! number of observed data at time t
                ENDIF
            ENDDO
c
c           if all data are missing at date t -> exit
            IF (mt .EQ. p) THEN
                info = -9001
                RETURN
            ENDIF
            IF ((mt .GT. 0).AND. (nt .GT. 0)) THEN ! if missing data
c
c           construct cov. matrix for miss and observed data
c           cov[missing,obs] : vector of missing vs. observed covariances
c           cov[obs,obs]     : matrix of observed covariances
c           cov[missing,missing] : matrix of missing covariances
            DO i = 1,mt
                k = iwork(pindMiss + i - 1)
                DO j = 1,nt
                    l = iwork(pindObs + j - 1)
                    dwork(pdcovMO + mt*(j-1) + i-1) =
     &              dwork(pdcov + p*(k-1) + l-1)
                ENDDO
            ENDDO
            CALL YMCPI ( p, dwork(pdcov), mt, iwork(pindMiss),
     &                   dwork(pdcovMM), info)
            CALL YMCPI ( p, dwork(pdcov), nt, iwork(pindObs),
     &                   dwork(pdcovOO), info)
            CALL JMS (nt,dwork(pdcovOO),dwork(pdjms),dwork(pdinvCovOO),
     &                info)
            IF (info .NE. 0) THEN
                info = -108 ! matrix is not sdp
                RETURN
            ENDIF
c
c           y(Miss) = M(Miss)+S(Miss,Obs)*inv(S(Obs,Obs))*(y(Obs)-M(Obs)
            CALL DV(nt,dwork(pdyObs),dwork(pdrhoObs),dwork(pdM))
            CALL PMV(nt,nt,dwork(pdinvCovOO),dwork(pdM),dwork(pdyM))
            CALL PMV(mt,nt,dwork(pdcovMO),dwork(pdyM),dwork(pdv))
            CALL SV(mt,dwork(pdrhoMiss),dwork(pdv),dwork(pdyMiss))

c
c           c(Miss,Miss) = S(Miss,Miss)-S(Miss,Obs)*inv(S(Obs,Obs))*S(Obs,Miss)
c           C(m,n) =  A(m,k) * B(k,n)
            CALL PM (mt,nt,nt,dwork(pdcovMO),dwork(pdinvCovOO),
     &                dwork(pdM1))
c           C(m,n) = A(m,k) * B'(n,k)
            CALL PMMT (mt,nt,mt,dwork(pdM1),dwork(pdcovMO),dwork(pdM2))
c           C(m,n) = A(m,n) - B(m,n)
            CALL DM (mt,mt,dwork(pdcovMM),dwork(pdM2),dwork(pdM3))
            DO i = 1,mt
                k = iwork(pindMiss + i - 1)
                DO j = 1,mt
                    l = iwork(pindMiss + j - 1)
                    dwork(pdCC + p*(k-1) + l - 1) =
     &              dwork(pdM3 + mt*(i-1) + j - 1)

                ENDDO
            ENDDO
            CALL SV (p*p, dwork(pdCC), dwork(pdMt), dwork(pdC))
            CALL YV (p*p, dwork(pdC), dwork(pdMt))
c
            mtt = 0
            DO i = 1,p
                IF (X(t,i) .LE. miss) THEN
                    mtt = mtt + 1
c                   replace missing data
                    dwork(pdy + mtt - 1) =
     &              MAX(dwork(pdyMiss + mtt - 1),0.0)
	              Y(t,i) = dwork(pdy + mtt - 1)
	            ENDIF
      		ENDDO
         ENDIF
       ENDDO
c
c       Step 2: update
c       update mean rho_new = mean(Y)
        CALL MCM (n,p,Y,dwork(pdrhoN))
        tmp = 1.0/n
c       update cov. matrix
        CALL IVX (p*p, dwork(pdM1), ZERO) ! M1 = zeros(p,p)
        CALL IVX (p*p, dwork(pdM2), ZERO) ! M2 = zeros(p,p)
        CALL IVX (p*p, dwork(pdM3), ZERO) ! M3 = zeros(p,p)
        CALL PVX (p*p, dwork(pdC), tmp, dwork(pdM1))
        CALL COVM (n,p,Y,dwork(pdcovm),dwork(pdM2),info)
        tmp = (n-1.0)/n
        CALL PVX (p*p, dwork(pdM2), tmp, dwork(pdM3))
        CALL SV (p*p, dwork(pdM1), dwork(pdM3), dwork(pdcovN))
c
c       test of algorithm convergence (mean and covariance matrix)
        CALL IVX (p, dwork(pdv), ZERO)    ! v = zeros(p)
        CALL IVX (p*p, dwork(pdM1), ZERO) ! M1 = zeros(p,p)
        CALL DV ( p, dwork(pdrho), dwork(pdrhoN), dwork(pdv) )
        CALL DV ( p*p, dwork(pdcov), dwork(pdcovN), dwork(pdM1) )
c       D4=[(M_new - M).^4; diag((S_new - S).^2) ]
        tmp = dwork(pdv) + dwork(pdM1)
        DO i = 2,p
            tmp = tmp + dwork(pdv+i-1)**4 + dwork(pdM1+p*(i-1)+i-1)**2
        ENDDO
c       Distance= mean(D4.^(1/4))
        Distance = SQRT(tmp/(2.0*p))
        Distance = SQRT(Distance)
        IF (Distance .LT. Tolerance) Convergence = 1
      ENDDO ! end while
c
      RETURN
      END

ext: f

max_stars_repo_name: alpgodev/numx

max_issues_repo_name: alpgodev/numx

max_forks_repo_name: alpgodev/numx


text:
Recently, a barrage of real people have been infected by the Front Page Davis Wiki Sockpuppet Virus. Were not quite sure how it is transferred from wiki user to wiki user, but we know one thing for sure: the icon is darn cute.

 Sockpuppet Virus Spreading Documentation

If you have been infected by the virus, document how you received it so we can find out who might be the Sockpuppet Mary.

  Sockpuppet Mary  No one knows their real identity.
   ?? Unknown Chain ??
    Users/JabberWokky  Self inflicted?
     Users/WesHardaker when Users/JabberWokky intentionally infected his Users/WesHardaker/Testing page.
      Users/EdHenn Ed Henns socks smell funny now because Users/WesHardaker had to pass it on to someone after all...
       Users/PhilipNeustrom Philip Wilt Chamberlain Neustrom contracted the virus from Users/EdHenn during one of his 10,000 Wiki encounters...  Then, he was miraculously cured!  Alas, he was left with a SockPuppet/Talk tragic physical condition that has adversely affected his ability to edit.
    Users/WilliamLewis  Very likely self inflicted.
    Users/DavidGrundler  Bought a supply to infect the below?  The tricky thing was how he affected Arlen after Arlen infected Jason.  We suspect a time machine was involved here.  (Or is it that the virus has a random amount of time where it is active but hasnt manifested itself.  Some can be infected for weeks without showing the telltale sign, while others exhibit the telltale sign immediately upon infection...)
     Users/ArlenAbraham
      Users/JasonAller
 Preventing Infection

As mentioned, how it spreads is unknown.  However, staying away from infected users will likely help.  Always practice safe editing!

 Origin

While the origin is unknown, it is worth noting that the virus appeared shortly after the Philip Neustrom Theory Neustromnomic theory Philip Neustrom Theory Number One number one was violated.  It is possible that the wiki itself is generating this darn cute icon to balance reality.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
 Prior Meetings (20052006)

  Wednesday Nov. 2nd at Cafe Roma by Campus at 4PM.
  Wednesday Nov. 9 at Cafe Roma by Campus at 5PM.
  Wednesday Nov. 16 at Cafe Roma by Campus at 5PM.
  Wednesday Nov. 23 at Cafe Roma by Campus at 4PM, photos available from http://www.broadbandreports.com/forum/remark,14871871 A. Leonard
  Wednesday Nov. 30 at Cafe Roma by Campus at 3:30PM.
  Wednesday Dec. 7 at Cafe Roma by Campus at 4PM.
  Wednesday Dec. 14 at Cafe Roma by Campus at 3:30PM.
  Wednesday Dec 21 at Cafe Roma by Campus at 4:00PM.
  Wednesday Dec 28 at Cafe Roma by Campus at 4:00PM.
  Wednesday January 4 at Cafe Roma by Campus at 3:30PM.
  Tuesday January 10 at Cafe Roma by Campus at 3:30PM.
  Tuesday January 17th at Cafe Roma by Campus, at 3:30PM.
  Saturday January 21st at Cafe Roma by Campus at 2:30PM. We took a drive to Sacramento where we photographed oldtown, the bridge, and the State Capitol.
  Tuesday January 24 at Cafe Roma by Campus at 3:30PM.
  Tuesday, February 7 at Cafe Roma by Campus at 3:30PM.
  Tuesday, February 14 at Cafe Roma by Campus at 3:30PM.
  Tuesday, February 21 at Cafe Roma by Campus at 3:30PM.
   Tuesday, February 28 at Cafe Roma by Campus at 3:30PM.
   Tuesday, March 7 at Cafe Roma by Campus at 4:30PM.
  Tuesday, March 14 at Cafe Roma by Campus at 3:15PM.


text:
      subroutine   iinit   ( n, a, x, incx )
C
C     ==================================================================
C
C     PURPOSE ... INITIALIZES INTEGER VECTOR TO
C                 A CONSTANT VALUE 'A'
C
C     CREATED ... AUG. 06, 1992
C
C     ==================================================================
C
C     -------------
C     ... ARGUMENTS
C     -------------
C
      INTEGER             N, A, X (*), INCX
C
C     -------------------
C     ... LOCAL VARIABLES
C     -------------------
C
      INTEGER             XADDR, I
C
C     ==================================================================
C
      IF  ( INCX .EQ. 1 )  THEN
C
C         ----------------------------------
C         ... UNIT INCREMENT (STANDARD CASE)
C         ----------------------------------
C
          DO 100 I = 1, N
              X(I) = A
  100     CONTINUE
C
      ELSE
C
C         ----------------------
C         ... NON-UNIT INCREMENT
C         ----------------------
C
          XADDR = 1
          IF  ( INCX .LT. 0 )  THEN
              XADDR = (-N+1)*INCX + 1
          ENDIF
C
          DO 200 I = 1, N
              X (XADDR) = A
              XADDR     = XADDR + INCX
  200     CONTINUE
C
      ENDIF
C
      RETURN
C
      END


text:
true constant dev

include ramen/ramen.f
ide

: platformer  s" ld sample/platformer/platformer" evaluate ;

permanent off  \ "gild" system assets
gild           \ and gild the dictionary

s" session.f" file-exists [if]
    include session.f
[then]

max_stars_repo_name: RogerLevy/RamenEngine

max_issues_repo_name: Alim-Oezdemir/RamenEngine

max_forks_repo_name: Alim-Oezdemir/RamenEngine


text:
program read_code
! 03/19/2022 09:45 AM allows the user to run executable again without recompiling
! 03/18/2022 10:25 PM read Fortran code from terminal, compile and run it
! to do:
! number the lines of code in the listing -- also have a "barelist"
! store the lines of code in an array of strings
! let user include code from another file
! let user delete a range of lines, including n1:n2, n: or :n
! let user edit a line
! let user insert a line
! let user comment out or uncomment a set of lines
! make it work on Linux
! allow the code to be linked with specified object or library files to create an executable
! indent code in if blocks and do loops
! expand p x,y to print*,x,y
! add end associate statements as needed
! check if each entered line of code is invalid and reject it if it is -- are there tools for this?
! let user make a block out of a group of lines
! allow tools such as Met Office stylist to analyze the code, or other tools to indent it
! move use statements and declarations with :: automatically to the write place
implicit none
character (len=*), parameter :: src = "temp.f90", full_src="junk.f90", &
                                executable = "a.exe"
character (len=1000)         :: text, code_line
integer, parameter           :: iu = 10, ju = 11, ku = 12
integer                      :: ierr,stat
open (unit=iu,file=full_src,iostat=stat)
if (stat == 0) close (iu,status="delete")
open (unit=iu,file=src,action="write",status="replace")
do
   write (*,"('> ')",advance="no")
   read (*,"(a)") text
   if (text == "quit" .or. text == "q") then
      stop "done -- code saved in " // trim(full_src)
   else if (text == "help" .or. text == "h") then
      print*,"quit or q to quit"
      print*,"run or r to compile and run executable"
      print*,"runold or ro to run executable again without recompiling"
      print*,"compile to compile without running executable"
      print*,"list or l for a program listing"
      print*,"clear or c to delete code and start over"
   else if (text == "clear" .or. text == "c") then
      close (iu,status="delete")
      open (unit=iu,file=src,action="write",status="replace")
   else if (text == "list" .or. text == "l") then
      close (iu)
      open (unit=ju,file=src,action="read")
      do
         read (ju,"(a)",iostat=ierr) code_line
         if (ierr == 0) then
            write (*,*) trim(code_line)
         else
            exit
         end if
      end do
      close (ju)
      open (unit=iu,file=src,action="write",position="append")
   else if (text == "run" .or. text == "r") then
      call compile()
      if (file_exists(executable)) call execute_command_line(executable)
   else if (text == "compile") then
      call compile()
   else if (text == "runold" .or. text == "ro") then
      if (file_exists(executable)) then
         call execute_command_line(executable)
      else
         print "(a,' does not exist')",trim(executable)
      end if
   else
      write (iu,"(a)") text ! add text to source file
!      print "(a)",trim(text)
   end if
end do
contains
!
subroutine compile()
integer :: stat, ku
call execute_command_line("copy " // src // " " // full_src // " > nul")
call execute_command_line("echo end >> " // full_src)
open (newunit=ku,file=executable,iostat=stat)
if (stat == 0) close(ku,status="delete") ! delete executable if it exists
call execute_command_line("gfortran " // trim(full_src))
end subroutine compile
!
function file_exists(file_name) result(exists)
character (len=*), intent(in) :: file_name
logical                       :: exists
inquire(file=file_name,exist=exists)
end function file_exists
!
end program read_code

max_stars_repo_name: Beliavsky/FortranShell

max_issues_repo_name: Beliavsky/FortranShell

max_forks_repo_name: Beliavsky/FortranShell


! simple example program
	program hello_world

	write (*,*) "Hello world."

	end program hello_world


      INTEGER N_HK, N_RT, N_TPPT
      REAL, CLDIMENSION(1,6) :: HK_RANGE
      REAL, DIMENSION(1,5) :: RT_RANGE
      REAL, DIMENSION(1,3) :: TPPT_RANGE
      REAL, DIMENSION(6,5,3) :: CEI_TABLE

      N_HK   = 6
      N_RT   = 5
      N_TPPT = 3

      HK_RANGE(1,1) = 7.0E0/5.0E0
      HK_RANGE(1,2) = 8.0E0/5.0E0
      HK_RANGE(1,3) = 9.0E0/5.0E0
      HK_RANGE(1,4) = 2.0E0
      HK_RANGE(1,5) = 1.1E1/5.0E0
      HK_RANGE(1,6) = 1.2E1/5.0E0

      RT_RANGE(1,1) = 8.0E2
      RT_RANGE(1,2) = 1.6E3
      RT_RANGE(1,3) = 3.2E3
      RT_RANGE(1,4) = 6.4E3
      RT_RANGE(1,5) = 1.28E4

      TPPT_RANGE(1,2) = 5.0E1
      TPPT_RANGE(1,3) = 1.0E2

      CEI_TABLE(1,1,2) = 2.872680733645279E-1
      CEI_TABLE(1,1,3) = 7.749716029388073E-1
      CEI_TABLE(1,2,2) = 2.039644535775086E-1
      CEI_TABLE(1,2,3) = 5.924318712770024E-1
      CEI_TABLE(1,3,2) = 1.443700895511554E-1
      CEI_TABLE(1,3,3) = 4.633605499494066E-1
      CEI_TABLE(1,4,2) = 1.004911816583118E-1
      CEI_TABLE(1,4,3) = 3.699178360235265E-1
      CEI_TABLE(1,5,2) = 6.733821510084832E-2
      CEI_TABLE(1,5,3) = 3.001150710959782E-1
      CEI_TABLE(2,1,2) = 2.15679590997346E-1
      CEI_TABLE(2,1,3) = 7.032901821333088E-1
      CEI_TABLE(2,2,2) = 1.38804898233508E-1
      CEI_TABLE(2,2,3) = 5.190092839824512E-1
      CEI_TABLE(2,3,2) = 8.540581665774229E-2
      CEI_TABLE(2,3,3) = 3.955784955811628E-1
      CEI_TABLE(2,4,2) = 4.687087515554092E-2
      CEI_TABLE(2,4,3) = 3.083346385006999E-1
      CEI_TABLE(2,5,2) = 1.823657536308939E-2
      CEI_TABLE(2,5,3) = 2.440957040061711E-1
      CEI_TABLE(3,1,2) = 1.351465804725244E-1
      CEI_TABLE(3,1,3) = 5.944539142560507E-1
      CEI_TABLE(3,2,2) = 6.948210676888233E-2
      CEI_TABLE(3,2,3) = 4.287398277454012E-1
      CEI_TABLE(3,3,2) = 2.446511411233687E-2
      CEI_TABLE(3,3,3) = 3.186820261578073E-1
      CEI_TABLE(3,4,3) = 2.413596674968999E-1
      CEI_TABLE(3,5,3) = 1.847296324018285E-1
      CEI_TABLE(4,1,2) = 5.497866179885413E-2
      CEI_TABLE(4,1,3) = 4.803653864888552E-1
      CEI_TABLE(4,2,2) = 1.021202577184958E-3
      CEI_TABLE(4,2,3) = 3.363622877434947E-1
      CEI_TABLE(4,3,3) = 2.409694933687837E-1
      CEI_TABLE(4,4,3) = 1.741794817845234E-1
      CEI_TABLE(4,5,3) = 1.254740804951762E-1
      CEI_TABLE(5,1,3) = 3.674090874850944E-1
      CEI_TABLE(5,2,3) = 2.450395028559258E-1
      CEI_TABLE(5,3,3) = 1.641449222712615E-1
      CEI_TABLE(5,4,3) = 1.077346391909527E-1
      CEI_TABLE(5,5,3) = 6.682986734578829E-2
      CEI_TABLE(6,1,3) = 2.566658710465144E-1
      CEI_TABLE(6,2,3) = 1.54999205732387E-1
      CEI_TABLE(6,3,3) = 8.805729195880023E-2
      CEI_TABLE(6,4,3) = 4.168825812144024E-2
      CEI_TABLE(6,5,3) = 8.364384007724771E-3


text:
module mod_test_functional_utilities
#include <messenger.h>

    use hdf5
    use h5lt

    use mod_hdf_utilities,         only: open_file_hdf, create_functional_group_hdf,                &
                                         create_functional_hdf, set_functional_reference_geom_hdf,  &
                                         set_functional_auxiliary_geom_hdf, set_functional_LS_hdf,  &
                                         close_file_hdf, remove_functional_hdf,                     &
                                         open_functional_group_hdf, close_functional_group_hdf

    implicit none

contains

    !>  Add a functional to a test mesh file.
    !!
    !!  @author Matteo Ugolotti
    !!  @date   12/12/2017
    !!
    !!  @param[in]  filename            Name of the mesh to which we need to add the functional
    !!  @param[in]  functional_name     The name of the registered functional to be added
    !!  @param[in]  reference_geom      Name of the reference geometry
    !!  @param[in]  auxiliary_geom      Name of the auxiliary geometry
    !!  @param[in]  linear_solver       Name of the linear solver
    !!
    !!
    !---------------------------------------------------------------------------------------------
    subroutine meshfile_add_functional(filename,functional_name,reference_geom,auxiliary_geom,linear_solver)
        character(*),                           intent(in)  :: filename
        character(*),                           intent(in)  :: functional_name
        character(*),                           intent(in)  :: reference_geom
        character(*),    optional,              intent(in)  :: auxiliary_geom
        character(*),    optional,              intent(in)  :: linear_solver

        integer(HID_T)     :: fid

        ! Open hdf meshfile
        fid = open_file_hdf(filename)

        ! Create the functional group if necessary
        call create_functional_group_hdf(fid)

        ! Create functional
        call create_functional_hdf(fid,functional_name)

        ! Define reference geometry
        call set_functional_reference_geom_hdf(fid,functional_name,reference_geom)

        ! Define auxiliary geometry
        if (present(auxiliary_geom)) call set_functional_auxiliary_geom_hdf(fid,functional_name,auxiliary_geom)

        ! Define linear solver
        if (present(linear_solver)) call set_functional_LS_hdf(fid,functional_name,linear_solver)

        ! Close hdf file
        call close_file_hdf(fid)

    end subroutine meshfile_add_functional
    !---------------------------------------------------------------------------------------------

    !>  Edit a functional to a test mesh file.
    !!
    !!  @author Matteo Ugolotti
    !!  @date   12/12/2017
    !!
    !!  @param[in]  filename            Name of the mesh to which we need to add the functional
    !!  @param[in]  functional_name     The name of the registered functional to be added
    !!  @param[in]  attribute           Attribute to be changed
    !!  @param[in]  new_entry           New entry
    !!
    !!
    !---------------------------------------------------------------------------------------------
    subroutine meshfile_edit_functional(filename,functional_name,attribute,new_entry)
        character(*),                           intent(in)  :: filename
        character(*),                           intent(in)  :: functional_name
        character(*),                           intent(in)  :: attribute
        character(*),                           intent(in)  :: new_entry

        integer(HID_T)     :: fid
        character(:), allocatable   :: user_msg

        ! Open hdf meshfile
        fid = open_file_hdf(filename)

        select case (trim(attribute))
            case("reference geometry")

                ! Define reference geometry
                call set_functional_reference_geom_hdf(fid,functional_name,new_entry)

            case("auxiliary geometry")

                ! Define auxiliary geometry
                call set_functional_auxiliary_geom_hdf(fid,functional_name,new_entry)

            case("linear solver")

                ! Define linear solver
                call set_functional_LS_hdf(fid,functional_name,new_entry)

            case default
                user_msg = "edit_functional2meshfile: There was no valid case that matched the incoming string"
                call chidg_signal(FATAL,user_msg)

        end select

        ! Close hdf file
        call close_file_hdf(fid)

    end subroutine meshfile_edit_functional
    !---------------------------------------------------------------------------------------------

    !>  Remove a functional from a test mesh file.
    !!
    !!  @author Matteo Ugolotti
    !!  @date   12/12/2017
    !!
    !!  @param[in]  filename            Name of the mesh to which we need to remove the functional
    !!  @param[in]  functional_name     The name of the registered functional to be removed
    !!
    !!
    !---------------------------------------------------------------------------------------------
    subroutine meshfile_remove_functional(filename,functional_name)
        character(*),                           intent(in)  :: filename
        character(*),                           intent(in)  :: functional_name

        integer(HID_T)     :: fid, fcl_id
        character(:), allocatable   :: user_msg

        ! Open hdf meshfile
        fid = open_file_hdf(filename)
        fcl_id = open_functional_group_hdf(fid)

        ! Remove functional
        call remove_functional_hdf(fcl_id,trim(functional_name))

        ! Close functional group
        call close_functional_group_hdf(fcl_id)
        ! Close hdf file
        call close_file_hdf(fid)

    end subroutine meshfile_remove_functional
    !---------------------------------------------------------------------------------------------

end module mod_test_functional_utilities

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG


text:
!-----------------------------------------------------------------------
!     Module:        diagno_bfield
!     Authors:       S. Lazerson (lazerson@pppl.gov)
!     Date:          03/02/2012
!     Description:   This subroutine calculates the field at a point in
!                    space.
!-----------------------------------------------------------------------
      SUBROUTINE diagno_bfield
!-----------------------------------------------------------------------
!     Libraries
!-----------------------------------------------------------------------
      USE stel_kinds, ONLY: rprec
      USE diagno_runtime
      USE virtual_casing_mod
      USE biotsavart
      USE safe_open_mod
      USE mpi_params
      USE mpi_inc
!-----------------------------------------------------------------------
!     Local Variables
!          ier            Error Flag
!          iunit          File ID Number
!-----------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER, PARAMETER :: BYTE_8 = SELECTED_INT_KIND (8)
#if defined(MPI_OPT)
      INTEGER(KIND=BYTE_8),ALLOCATABLE :: mnum(:), moffsets(:)
      INTEGER :: numprocs_local, mylocalid, mylocalmaster
      INTEGER :: MPI_COMM_LOCAL
#endif
      INTEGER(KIND=BYTE_8) :: icount, chunk
      INTEGER :: ier, iunit, ncoils, i, ig, iunit_out, ncg
      REAL(rprec), ALLOCATABLE, DIMENSION(:) :: xp, yp, rp, phip, zp, bx, by, br, bphi, bz, modb,&
                     bxp, byp, bzp
      REAL(rprec), ALLOCATABLE, DIMENSION(:,:) :: bx_mut, by_mut, bz_mut
      REAL(rprec) :: xvec(3),bvec(3)

!-----------------------------------------------------------------------
!     Begin Subroutine
!-----------------------------------------------------------------------
      ! Basic copy of MPI_COMM_DIANGO
      mylocalid = myworkid
      mylocalmaster = master
      MPI_COMM_LOCAL = MPI_COMM_DIAGNO
      numprocs_local = nprocs_diagno
      ncg = SIZE(coil_group)

      ! Read File
      if(lverb) write(6,*)' ---Calculating Magnetic Field Test Points'
      IF(lverb .and. lrphiz) THEN
         WRITE(6,'(14X,A,7(6X,A,6X))') 'No','r','phi','z','Br','Bphi','Bz','|B|'
      ELSE IF (lverb) THEN
         WRITE(6,'(14X,A,7(6X,A,6X))') 'No','x','y','z','Bx','By','Bz','|B|'
      END IF
      iunit = 26; iunit_out = 27;
      IF (myworkid == master) THEN
         CALL safe_open(iunit,ier,TRIM(bfield_points_file),'old','formatted')
         READ(iunit,*) ncoils
         ALLOCATE(xp(ncoils), yp(ncoils), rp(ncoils), phip(ncoils), zp(ncoils), bx(ncoils), &
                  by(ncoils), br(ncoils), bphi(ncoils), bz(ncoils), modb(ncoils),&
                  bxp(ncoils), byp(ncoils), bzp(ncoils))
         DO i = 1, ncoils
            IF (lrphiz) THEN
               READ(iunit,*) rp(i),phip(i),zp(i)
               xp(i) = rp(i) * cos(phip(i))
               yp(i) = rp(i) * sin(phip(i))
            ELSE
               READ(iunit,*) xp(i),yp(i),zp(i)
            END IF
         END DO
         CLOSE(iunit)
      END IF

#if defined(MPI_OPT)
      CALL MPI_BARRIER(MPI_COMM_DIAGNO,ierr_mpi)
      IF (ierr_mpi /=0) CALL handle_err(MPI_BARRIER_ERR,'diagno_bfield1',ierr_mpi)
      CALL MPI_BCAST(ncoils,1,MPI_INTEGER, master, MPI_COMM_DIAGNO,ierr_mpi)
      IF (ierr_mpi /=0) CALL handle_err(MPI_BCAST_ERR,'diagno_bfield1',ierr_mpi)
#endif
      IF (myworkid /= master) ALLOCATE(xp(ncoils), yp(ncoils), rp(ncoils), phip(ncoils), zp(ncoils), bx(ncoils), &
                                       by(ncoils), br(ncoils), bphi(ncoils), bz(ncoils), modb(ncoils),&
                                       bxp(ncoils), byp(ncoils), bzp(ncoils))
      ALLOCATE(bx_mut(ncoils,ncg), by_mut(ncoils,ncg), bz_mut(ncoils,ncg))

#if defined(MPI_OPT)
      CALL MPI_BARRIER(MPI_COMM_DIAGNO,ierr_mpi)
      IF (ierr_mpi /=0) CALL handle_err(MPI_BARRIER_ERR,'diagno_bfield2',ierr_mpi)
      CALL MPI_BCAST(xp,ncoils,MPI_DOUBLE_PRECISION, master, MPI_COMM_DIAGNO,ierr_mpi)
      IF (ierr_mpi /=0) CALL handle_err(MPI_BCAST_ERR,'diagno_bfield2',ierr_mpi)
      CALL MPI_BCAST(yp,ncoils,MPI_DOUBLE_PRECISION, master, MPI_COMM_DIAGNO,ierr_mpi)
      IF (ierr_mpi /=0) CALL handle_err(MPI_BCAST_ERR,'diagno_bfield2',ierr_mpi)
      CALL MPI_BCAST(zp,ncoils,MPI_DOUBLE_PRECISION, master, MPI_COMM_DIAGNO,ierr_mpi)
      IF (ierr_mpi /=0) CALL handle_err(MPI_BCAST_ERR,'diagno_bfield2',ierr_mpi)
      CALL MPI_BCAST(rp,ncoils,MPI_DOUBLE_PRECISION, master, MPI_COMM_DIAGNO,ierr_mpi)
      IF (ierr_mpi /=0) CALL handle_err(MPI_BCAST_ERR,'diagno_bfield2',ierr_mpi)
      CALL MPI_BCAST(phip,ncoils,MPI_DOUBLE_PRECISION, master, MPI_COMM_DIAGNO,ierr_mpi)
      IF (ierr_mpi /=0) CALL handle_err(MPI_BCAST_ERR,'diagno_bfield2',ierr_mpi)
#endif

      ! Divide up the work
      chunk = FLOOR(REAL(ncoils) / REAL(numprocs_local))
      mystart = myworkid*chunk + 1
      myend = mystart + chunk - 1
#if defined(MPI_OPT)
      IF (ALLOCATED(mnum)) DEALLOCATE(mnum)
      IF (ALLOCATED(moffsets)) DEALLOCATE(moffsets)
      ALLOCATE(mnum(numprocs_local), moffsets(numprocs_local))
      CALL MPI_ALLGATHER(chunk,1,MPI_INTEGER,mnum,1,MPI_INTEGER,MPI_COMM_LOCAL,ierr_mpi)
      CALL MPI_ALLGATHER(mystart,1,MPI_INTEGER,moffsets,1,MPI_INTEGER,MPI_COMM_LOCAL,ierr_mpi)
      i = 1
      DO
         IF ((moffsets(numprocs_local)+mnum(numprocs_local)-1) == ncoils) EXIT
         IF (i == numprocs_local) i = 1
         mnum(i) = mnum(i) + 1
         moffsets(i+1:numprocs_local) = moffsets(i+1:numprocs_local) + 1
         i=i+1
      END DO
      mystart = moffsets(mylocalid+1)
      chunk  = mnum(mylocalid+1)
      myend   = mystart + chunk - 1
#endif

      ! Do Work
      bx = 0; by = 0; bz = 0;
      bx_mut = 0; by_mut = 0; bz_mut = 0;
      DO i = mystart, myend
         xvec(1) = xp(i); xvec(2)=yp(i); xvec(3)=zp(i)
         IF (lcoil) THEN
            DO ig = 1, ncg
               bvec=0
               CALL bsc_b(coil_group(ig),xvec,bvec)
               bx_mut(i,ig) = bvec(1)
               by_mut(i,ig) = bvec(2)
               bz_mut(i,ig) = bvec(3)
            END DO
         END IF
         IF (lmut .or. lvac) CYCLE
         CALL bfield_vc(xp(i),yp(i),zp(i),bx(i),by(i),bz(i),ier)
      END DO

      ! Now handle the arrays
      IF (lmut) THEN
         IF (myworkid == master) THEN
            CALL MPI_REDUCE(MPI_IN_PLACE,bx_mut,ncoils*ncg,MPI_DOUBLE_PRECISION,MPI_SUM,master,MPI_COMM_DIAGNO,ierr_mpi)
            CALL MPI_REDUCE(MPI_IN_PLACE,by_mut,ncoils*ncg,MPI_DOUBLE_PRECISION,MPI_SUM,master,MPI_COMM_DIAGNO,ierr_mpi)
            CALL MPI_REDUCE(MPI_IN_PLACE,bz_mut,ncoils*ncg,MPI_DOUBLE_PRECISION,MPI_SUM,master,MPI_COMM_DIAGNO,ierr_mpi)
         ELSE
            CALL MPI_REDUCE(bx_mut,bx_mut,ncoils*ncg,MPI_DOUBLE_PRECISION,MPI_SUM,master,MPI_COMM_DIAGNO,ierr_mpi)
            CALL MPI_REDUCE(by_mut,by_mut,ncoils*ncg,MPI_DOUBLE_PRECISION,MPI_SUM,master,MPI_COMM_DIAGNO,ierr_mpi)
            CALL MPI_REDUCE(bz_mut,bz_mut,ncoils*ncg,MPI_DOUBLE_PRECISION,MPI_SUM,master,MPI_COMM_DIAGNO,ierr_mpi)
         END IF
      ELSE
         ! Plasma + vacuum
         bx = bx + SUM(bx_mut, DIM=2)
         by = by + SUM(by_mut, DIM=2)
         bz = bz + SUM(bz_mut, DIM=2)
#if defined(MPI_OPT)
         CALL MPI_BARRIER(MPI_COMM_LOCAL,ierr_mpi)
         IF (ierr_mpi /=0) CALL handle_err(MPI_BARRIER_ERR,'diagno_bfield3',ierr_mpi)
         CALL MPI_ALLGATHERV(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,&
                        bx,mnum,moffsets-1,MPI_DOUBLE_PRECISION,&
                        MPI_COMM_LOCAL,ierr_mpi)
         CALL MPI_ALLGATHERV(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,&
                        by,mnum,moffsets-1,MPI_DOUBLE_PRECISION,&
                        MPI_COMM_LOCAL,ierr_mpi)
         CALL MPI_ALLGATHERV(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,&
                        bz,mnum,moffsets-1,MPI_DOUBLE_PRECISION,&
                        MPI_COMM_LOCAL,ierr_mpi)
#endif
      END IF

      ! PRINT TO SCREEN
      IF (myworkid == master) THEN
         IF (lmut) THEN
         ELSE
            CALL safe_open(iunit_out,ier,'diagno_btest.'//TRIM(id_string),'replace','formatted')
            modb = sqrt(bx*bx+by*by+bz*bz)
            br   = bx * cos(phip) + by * sin(phip)
            bphi = by * cos(phip) - bx * sin(phip)
            IF (lrphiz) THEN
               DO i = 1, ncoils
                  IF (lverb) WRITE(6,'(13X,I8,1X,7E14.5)') i,rp(i),phip(i),zp(i),br(i),bphi(i),bz(i),modb(i)
                  WRITE(iunit_out,'(13X,I8,1X,7E14.5)') i,rp(i),phip(i),zp(i),br(i),bphi(i),bz(i),modb(i)
               END DO
               WRITE(iunit_out,'(A)')'   #   rp[m]    phip[deg]    zp[m]      B_R[T]    B_PHI[T]      B_Z[T]      |B|[T]'
            ELSE
               DO i = 1, ncoils
                  IF (lverb) WRITE(6,'(13X,I8,1X,7E14.5)') i,xp(i),yp(i),zp(i),bx(i),by(i),bz(i),modb(i)
                  WRITE(iunit_out,'(13X,I8,1X,7E14.5)') i,xp(i),yp(i),zp(i),bx(i),by(i),bz(i),modb(i)
               END DO
               WRITE(iunit_out,'(A)')'   #   xp[m]      yp[m]      zp[m]      B_X[T]      B_Y[T]      B_Z[T]      |B|[T]'
            END IF
            CLOSE(iunit_out)
         END IF
      END IF

      ! Clean up
      DEALLOCATE(xp, yp, rp, phip, zp, bx, by, br, bphi, bz, modb,bxp, byp, bzp)
      DEALLOCATE(bx_mut,by_mut,bz_mut)
      RETURN
!-----------------------------------------------------------------------
!     End Subroutine
!-----------------------------------------------------------------------
      END SUBROUTINE diagno_bfield

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
!------------------------------------------------------------------------!
!  The Community Multiscale Air Quality (CMAQ) system software is in     !
!  continuous development by various groups and is based on information  !
!  from these groups: Federal Government employees, contractors working  !
!  within a United States Government contract, and non-Federal sources   !
!  including research institutions.  These groups give the Government    !
!  permission to use, prepare derivative works of, and distribute copies !
!  of their work in the CMAQ system to the public and to permit others   !
!  to do so.  The United States Environmental Protection Agency          !
!  therefore grants similar permission to use the CMAQ system software,  !
!  but users are requested to provide copies of derivative works or      !
!  products designed to operate in the CMAQ system to the United States  !
!  Government without restrictions as to use by others.  Software        !
!  that is used with the CMAQ system but distributed under the GNU       !
!  General Public License or the GNU Lesser General Public License is    !
!  subject to their copyright restrictions.                              !
!------------------------------------------------------------------------!

C RCS file, release, date & time of last delta, author, state, [and locker]
C $Header: /project/work/rep/STENEX/src/se_snl/se_domain_info_ext.f,v 1.1 2004/03/26 16:16:47 yoj Exp $

C what(1) key, module and SID; SCCS file; date and time of last delta:
C %W% %P% %G% %U%

C --------------------------------------------------------------------------
C Note: all these variables with prefix se_ are for stencil exchange library
C       only
C
C to define domain decomposition info variables:
C
C -- global values
C
C   se_gl_ncols  -- number of columns in the physical domain
C   se_gl_nrows  -- number of rows in the physical domain
C   se_gl_nlays  -- number of layers in the physical domain
C   se_gl_nspcs  -- number of species in the physical domain
C   se_gl_ind    -- global fine domain to processor map: the first two entries
C                   are the beginning and ending row number, the second two
C                   entries are the beginning and ending column number, second
C                   dimension corresponding to the processor number
C
C -- local dimension values
C
C   se_ncols -- max. number of columns in a processor for data declaration
C   se_nrows -- max. number of rows in a processor for data declaration
C   se_nlays -- max. number of layers in a processor for data declaration
C   se_nspcs -- max. number of species in a processor for data declaration
C
C -- local processor data dimension values
C
C   se_my_ncols -- number of columns a processor is responsible for
C   se_my_nrows -- number of rows a processor is responsible for
C   se_my_nlays -- number of layers a processor is responsible for
C   se_my_nspcs -- number of species a processor is responsible for
C
C -- local processor data dimension values plus one
C
C   se_my_ncolsp1 -- number of columns a processor is responsible for plus one
C   se_my_nrowsp1 -- number of rows a processor is responsible for plus one
C   se_my_nlaysp1 -- number of layers a processor is responsible for plus one
C   se_my_nspcsp1 -- number of species a processor is responsible for plus one
C
C eg. given eight processors, four processors are assigned to column
C     dimension and 2 are assigned to row dimension (note: vertical axis
C     denotes the column dimension and horizontal axis denotes row dimension),
C     and a 2-D domain of size 17 by 11
C
C     ie. numprocs = 8               PE map:    0   1
C         npcol    = 4                          2   3
C         nprow    = 2                          4   5
C         gl_ncols = 17                         6   7
C         gl_nrows = 11
C
C     in each PE, ncols by nrows = 5 x 6
C     se_my_ncols by se_my_nrows = 5 x 6 for PE 0
C     se_my_ncols by se_my_nrows = 5 x 5 for PE 1
C     se_my_ncols by se_my_nrows = 4 x 6 for PE 2, 4, and 6
C     se_my_ncols by se_my_nrows = 4 x 5 for PE 3, 5, and 7
C --------------------------------------------------------------------------

        module se_domain_info_ext

          integer :: se_gl_ncols
          integer :: se_gl_nrows
          integer :: se_gl_nlays
          integer :: se_gl_nspcs

          integer, allocatable, save, target :: se_gl_ind(:, :, :)
          integer, pointer :: se_gl_ind_ptr (:, :, :)

          integer :: se_ncols
          integer :: se_nrows
          integer :: se_nlays
          integer :: se_nspcs

          integer :: se_my_ncols
          integer :: se_my_nrows
          integer :: se_my_nlays
          integer :: se_my_nspcs
          integer :: se_my_ncolsp1
          integer :: se_my_nrowsp1
          integer :: se_my_nlaysp1
          integer :: se_my_nspcsp1

        end module se_domain_info_ext

max_stars_repo_name: Simeng-unique/CMAQ-changed

max_issues_repo_name: Simeng-unique/CMAQ-changed

max_forks_repo_name: Simeng-unique/CMAQ-changed


text:
subroutine opt_reanut(a,itask)
   !Read Numerical Treatment
   use typre
   use Mod_Listen
   use Mod_Mesh
   use Mod_Optics
   implicit none
   class (OpticsProblem) :: a
   integer(ip) :: itask

   integer(ip) :: istab,ipart,rpalo

   !For Listener
   real(rp), pointer     :: param(:) => NULL()
   character(5), pointer :: words(:) => NULL()
   integer(ip), pointer  :: nnpar => NULL(),nnwor => NULL()

   call a%Listener%getarrs(words,param,nnpar,nnwor)

   !Initializations (defaults)
   if (itask == 0) then

      a%kfl_dissi = 0   !Default dissipation scheme is Smagorinsky
      a%a_obukhov = 2.4

   !Inside the Numerical Treatment Block
   elseif (itask == 1) then
      if(words(1)=='DISSI') then
         if (a%Listener%exists('SMAGO')) then
            a%kfl_dissi = 0
            a%turbu(1)=a%Listener%param(2)
         elseif (a%Listener%exists('WALE ')) then
            a%kfl_dissi = 2
            a%turbu(1)=a%Listener%param(2)

         elseif (a%Listener%exists('EXTER')) then
            a%kfl_dissi = 1
         endif
      elseif(words(1)=='OBUKH') then
         a%a_Obukhov = param(1)
      endif

   !Finalize reading operations
   elseif (itask == 100) then

   endif

end subroutine opt_reanut

max_stars_repo_name: ciaid-colombia/InsFEM

max_issues_repo_name: ciaid-colombia/InsFEM

max_forks_repo_name: ciaid-colombia/InsFEM


text:
C$Procedure      TISBOD ( Transformation, inertial state to bodyfixed )

      SUBROUTINE TISBOD ( REF, BODY, ET, TSIPM )

C$ Abstract
C
C     Return a 6x6 matrix that transforms states in inertial
C     coordinates to states in body-equator-and-prime-meridian
C     coordinates.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     PCK
C     NAIF_IDS
C     ROTATION
C     TIME
C
C$ Keywords
C
C     TRANSFORMATION
C     ROTATION
C
C$ Declarations

      IMPLICIT NONE

      INCLUDE               'errhnd.inc'
      INCLUDE               'frmtyp.inc'
      INCLUDE               'zzctr.inc'

      CHARACTER*(*)         REF
      INTEGER               BODY
      DOUBLE PRECISION      ET
      DOUBLE PRECISION      TSIPM   ( 6,6 )

C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  --------------------------------------------------
C     REF        I   ID of inertial reference frame to transform from
C     BODY       I   ID code of body
C     ET         I   Epoch of transformation
C     TSIPM      O   Transformation (state), inertial to prime meridian
C
C$ Detailed_Input
C
C     REF         is the NAIF name for an inertial reference frame.
C                 Acceptable names include:
C
C                    Name       Description
C                    --------   --------------------------------
C                    'J2000'    Earth mean equator, dynamical
C                               equinox of J2000
C
C                    'B1950'    Earth mean equator, dynamical
C                               equinox of B1950
C
C                    'FK4'      Fundamental Catalog (4)
C
C                    'DE-118'   JPL Developmental Ephemeris (118)
C
C                    'DE-96'    JPL Developmental Ephemeris ( 96)
C
C                    'DE-102'   JPL Developmental Ephemeris (102)
C
C                    'DE-108'   JPL Developmental Ephemeris (108)
C
C                    'DE-111'   JPL Developmental Ephemeris (111)
C
C                    'DE-114'   JPL Developmental Ephemeris (114)
C
C                    'DE-122'   JPL Developmental Ephemeris (122)
C
C                    'DE-125'   JPL Developmental Ephemeris (125)
C
C                    'DE-130'   JPL Developmental Ephemeris (130)
C
C                    'GALACTIC' Galactic System II
C
C                    'DE-200'   JPL Developmental Ephemeris (200)
C
C                    'DE-202'   JPL Developmental Ephemeris (202)
C
C                 (See the routine CHGIRF for a full list of names.)
C
C                 The output TIPM will give the transformation
C                 from this frame to the bodyfixed frame specified by
C                 BODY at the epoch specified by ET.
C
C     BODY        is the integer ID code of the body for which the
C                 state transformation matrix is requested. Bodies
C                 are numbered according to the standard NAIF
C                 numbering scheme. The numbering scheme is
C                 explained in the NAIF_IDS required reading file.
C
C     ET          is the epoch at which the state transformation
C                 matrix is requested. (This is typically the
C                 epoch of observation minus the one-way light time
C                 from the observer to the body at the epoch of
C                 observation.)
C
C$ Detailed_Output
C
C     TSIPM       is a 6x6 transformation matrix. It is used to
C                 transform states from inertial coordinates to body
C                 fixed (also called equator and prime meridian ---
C                 PM) coordinates.
C
C                 Given a state S in the inertial reference frame
C                 specified by REF, the corresponding bodyfixed state
C                 is given by the matrix vector product:
C
C                    TSIPM * S
C
C                 The X axis of the PM system is directed  to the
C                 intersection of the equator and prime meridian.
C                 The Z axis points along  the spin axis and points
C                 towards the same side of the invariable plane of
C                 the solar system as does earth's north pole.
C
C                 NOTE: The inverse of TSIPM is NOT its transpose.
C                       The matrix, TSIPM, has a structure as shown
C                       below:
C
C                             -            -
C                            |       :      |
C                            |   R   :  0   |
C                            | ......:......|
C                            |       :      |
C                            | dR_dt :  R   |
C                            |       :      |
C                             -            -
C
C                       where R is a time varying rotation matrix and
C                       dR_dt is its derivative. The inverse of this
C                       matrix is:
C
C                             -              -
C                            |     T  :       |
C                            |    R   :  0    |
C                            | .......:.......|
C                            |        :       |
C                            |      T :   T   |
C                            | dR_dt  :  R    |
C                            |        :       |
C                             -              -
C
C                       The SPICE routine INVSTM is available for
C                       producing this inverse.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C     1) If data required to define the body-fixed frame associated
C        with BODY are not found in the binary PCK system or the kernel
C        pool, the error SPICE(FRAMEDATANOTFOUND) is signaled. In
C        the case of IAU style body-fixed frames, the absence of
C        prime meridian polynomial data (which are required) is used
C        as an indicator of missing data.
C
C     2) If the test for exception (1) passes, but in fact requested
C        data are not available in the kernel pool, the error will be
C        signaled by routines in the call tree of this routine.
C
C     3) If the kernel pool does not contain all of the data required
C        to define the number of nutation precession angles
C        corresponding to the available nutation precession
C        coefficients, the error SPICE(INSUFFICIENTANGLES) is
C        signaled.
C
C     4) If the reference frame REF is not recognized, a routine
C        called by TISBOD will diagnose the condition and invoke the
C        SPICE error handling system.
C
C     5) If the specified body code BODY is not recognized, the
C        error is diagnosed by a routine called by TISBOD.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     The matrix for transforming inertial states to bodyfixed
C     states is the 6x6 matrix shown below as a block structured
C     matrix.
C
C                 -            -
C                |       :      |
C                | TIPM  :  0   |
C                | ......:......|
C                |       :      |
C                | DTIPM : TIPM |
C                |       :      |
C                 -            -
C
C     This can also be expressed in terms of Euler angles
C     PHI, DELTA and W. The transformation from inertial to
C     bodyfixed coordinates is represented in the SPICE kernel
C     pool as:
C
C            TIPM = [W] [DELTA] [PHI]
C                      3       1     3
C     Thus
C
C           DTIPM = d[W] /dt [DELTA] [PHI]
C                       3           1     3
C
C                 + [W] d[DELTA] /dt [PHI]
C                      3             1     3
C
C                 + [W] [DELTA] d[PHI] /dt
C                      3       1          3
C
C
C     If a binary PCK file record can be used for the time and
C     body requested, it will be used. The most recently loaded
C     binary PCK file has first priority, followed by previously
C     loaded binary PCK files in backward time order. If no
C     binary PCK file has been loaded, the text P_constants
C     kernel file is used.
C
C     If there is only text PCK kernel information, it is
C     expressed in terms of RA, DEC and W (same W as above), where
C
C        RA    = PHI - HALFPI()
C        DEC   = HALFPI() - DELTA
C
C     The angles RA, DEC, and W are defined as follows in the
C     text PCK file:
C
C                                         2      ____
C                                    RA2*t       \
C            RA  = RA0  + RA1*t/T  + ------   +  /     a  sin theta
C                                       2        ----   i          i
C                                      T           i
C
C                                          2     ____
C                                    DEC2*t      \
C            DEC = DEC0 + DEC1*t/T + -------  +  /    d  cos theta
C                                        2       ----  i          i
C                                       T          i
C
C
C                                        2      ____
C                                    W2*t       \
C            W   = W0   + W1*t/d   + -----   +  /     w  sin theta
C                                       2       ----   i          i
C                                      d          i
C
C
C     where:
C
C            d = seconds/day
C
C            T = seconds/Julian century
C
C            a , d , and w  arrays apply to satellites only.
C             i   i       i
C
C            theta  = THETA0(i) + THETA1(i)*t/T are specific to each
C                 i
C
C            planet.
C
C
C     These angles -- typically nodal rates -- vary in number and
C     definition from one planetary system to the next.
C
C     Thus
C                                         ____
C                             2*RA2*t     \
C          dRA/dt  = RA1/T  + -------   + /   a THETA1(i)/T cos theta
C                                 2       ---- i                     i
C                                T          i
C
C                                         ____
C                             2*DEC2*t    \
C          dDEC/dt = DEC1/T + --------  - /    d  THETA1(i)/T sin theta
C                                 2       ----  i                      i
C                                T          i
C
C                                         ____
C                             2*W2*t      \
C          dW/dt   = W1/d   + ------    + /    w  THETA1(i)/T cos theta
C                                 2       ----  i                      i
C                                d          i
C
C
C$ Examples
C
C     Note that the data needed to compute the output state transition
C     matrix must have been made available to your program by having
C     loaded an appropriate binary or text PCK file via FURNSH.
C
C     Example 1.
C
C     In the following code fragment, TISBOD is used to transform
C     a state in J2000 inertial coordinates to a state in bodyfixed
C     coordinates.
C
C     The 6-vectors EULANG represents the inertial state (position and
C     velocity) of an object with respect to the center of the body
C     at time ET.
C
C     C
C     C     First load the kernel pool.
C     C
C           CALL FURNSH ( 'PLANETARY_CONSTANTS.KER' )
C
C     C
C     C     Next get the transformation and its derivative.
C     C
C           CALL TISBOD ( 'J2000', BODY, ET, TSIPM )
C
C     C
C     C     Convert position to bodyfixed coordinates.
C     C
C           CALL MXVG    ( TSIPM, EULANG, 6, 6, BDSTAT )
C
C
C     Example 2.
C
C     In the example below, TISBOD is used to compute
C     the angular velocity vector (with respect to an inertial frame)
C     of the specified body at time ET.
C
C     C
C     C     First get the state transformation matrix.
C     C
C           CALL TISBOD ( BODY,  ET,   TSIPM )
C
C     C
C     C     This matrix has the form:
C     C
C     C           -            -
C     C          |       :      |
C     C          | TIPM  :  0   |
C     C          | ......:......|
C     C          |       :      |
C     C          | DTIPM : TIPM |
C     C          |       :      |
C     C           -            -
C     C
C     C     We extract TIPM and DTIPM
C     C
C
C           DO  I = 1,3
C              DO  J = 1,3
C
C                 TIPM  ( I, J ) = TSIPM ( I,   J )
C                 DTIPM ( I, J ) = TSIPM ( I+3, J )
C
C              END DO
C           END DO
C
C     C
C     C     The transpose of TIPM and DTIPM, (TPMI and DTPMI), give
C     C     the transformation from bodyfixed coordinates to inertial
C     C     coordinates.
C     C
C     C     Here is a fact about the relationship between angular
C     C     velocity associated with a time varying rotation matrix
C     C     that gives the orientation of a body with respect to
C     C     an inertial frame.
C     C
C     C        The angular velocity vector can be read from the off
C     C        diagonal components of the matrix product:
C     C
C     C                                t
C     C        OMEGA =     DTPMI * TPMI
C     C
C     C                         t
C     C              =     DTIPM * TIPM
C     C
C     C        the components of the angular velocity V will appear
C     C        in this matrix as:
C     C
C     C             _                   _
C     C            |                     |
C     C            |   0    -V(3)  V(2)  |
C     C            |                     |
C     C            |  V(3)    0   -V(1)  |
C     C            |                     |
C     C            | -V(2)   V(1)   0    |
C     C            |_                   _|
C     C
C     C
C           CALL MTXM ( DTIPM, TIPM, OMEGA )
C
C           V(1) = OMEGA (3,2)
C           V(2) = OMEGA (1,3)
C           V(3) = OMEGA (2,1)
C
C$ Restrictions
C
C     The kernel pool must be loaded with the appropriate coefficients
C     (from the P_constants kernel or binary PCK file) prior to
C     calling this routine.
C
C$ Literature_References
C
C     None.
C
C$ Author_and_Institution
C
C     N. J. Bachman   (JPL)
C     B. V. Semenov   (JPL)
C     W. L. Taber     (JPL)
C     K. S. Zukor     (JPL)
C
C$ Version
C
C-    SPICELIB Version 4.5.0, 26-JUL-2016 (BVS)
C
C        The routine was updated to be more efficient by using a hash
C        and buffers so save text PCK data instead of doing kernel POOL
C        look-ups over an over again. The routine now checks the POOL
C        state counter and dumps all buffered data if it changes.
C
C        BUG FIX: changed available room in the BODVCD call
C        fetching 'NUT_PREC_ANGLES' from MAXANG to MAXANG*2.
C
C-    SPICELIB Version 4.4.0, 01-FEB-2008 (NJB)
C
C        The routine was updated to improve the error messages created
C        when required PCK data are not found. Now in most cases the
C        messages are created locally rather than by the kernel pool
C        access routines. In particular missing binary PCK data will
C        be indicated with a reasonable error message.
C
C-    SPICELIB Version 4.3.0, 13-DEC-2005 (NJB)
C
C        Bug fix:  previous update introduced bug in state
C        transformation when REF was unequal to PCK native frame.
C
C-    SPICELIB Version 4.2.0, 23-OCT-2005 (NJB)
C
C        Re-wrote portions of algorithm to simplify source code.
C        Updated to remove non-standard use of duplicate arguments
C        in MXM and VADDG calls.
C
C        Replaced calls to ZZBODVCD with calls to BODVCD.
C
C-    SPICELIB Version 4.1.0, 05-JAN-2005 (NJB)
C
C        Tests of routine FAILED() were added.
C
C-     SPICELIB Version 4.0.0, 12-FEB-2004 (NJB)
C
C         Code has been updated to support satellite ID codes in the
C         range 10000 to 99999 and to allow nutation precession angles
C         to be associated with any object.
C
C         Implementation changes were made to improve robustness
C         of the code.
C
C-     SPICELIB Version 3.3.0, 29-MAR-1995 (WLT)
C
C         Properly initialized the variable NPAIRS.
C
C-     SPICELIB Version 3.2.0, 22-MAR-1995 (KSZ)
C
C         Changed to call PCKMAT rather than PCKEUL.
C
C-     SPICELIB Version 3.1.0, 18-OCT-1994 (KSZ)
C
C         Fixed bug which incorrectly modded DW by two pi.
C
C-     SPICELIB Version 3.0.0, 10-MAR-1994 (KSZ)
C
C         Changed to look for binary PCK file, and used this
C         to find Euler angles, if such data has been loaded.
C
C-     SPICELIB Version 2.0.1, 10-MAR-1992 (WLT)
C
C         Comment section for permuted index source lines was added
C         following the header.
C
C-     SPICELIB Version 2.0.0, 04-SEP-1991 (NJB)
C
C         Updated to handle P_constants referenced to different epochs
C         and inertial reference frames.
C
C         $Required_Reading and $Literature_References sections were
C         updated.
C
C-     SPICELIB Version 1.0.0, 05-NOV-1990 (WLT)
C
C-&

C$ Index_Entries
C
C     transformation from inertial state to bodyfixed
C
C-&

C$ Revisions
C
C-    SPICELIB Version 4.5.0, 02-MAR-2016 (BVS)
C
C        The routine was updated to be more efficient by using a hash
C        and buffers so save text PCK data instead of doing kernel POOL
C        look-ups over an over again. The routine now checks the POOL
C        state counter and dumps all buffered data if it changes.
C
C-    SPICELIB Version 4.2.0, 06-SEP-2005 (NJB)
C
C        Re-wrote portions of algorithm to simplify source code.
C        The routine now takes advantage of EUL2XF, which wasn't
C        available when the first version of this routine was written.
C
C        Updated to remove non-standard use of duplicate arguments
C        in MXM and VADDG calls.
C
C        Replaced calls to ZZBODVCD with calls to BODVCD.
C
C-    SPICELIB Version 1.1.0, 05-JAN-2005 (NJB)
C
C        Tests of routine FAILED() were added. The new checks
C        are intended to prevent arithmetic operations from
C        being performed with uninitialized or invalid data.
C
C-     SPICELIB Version 4.0.0, 27-JAN-2004 (NJB)
C
C         Code has been updated to support satellite ID codes in the
C         range 10000 to 99999 and to allow nutation precession angles
C         to be associated with any object.
C
C         Calls to deprecated kernel pool access routine RTPOOL
C         were replaced by calls to GDPOOL.
C
C         Calls to BODVAR have been replaced with calls to
C         ZZBODVCD.
C
C-     SPICELIB Version 3.3.0, 29-MAR-1995 (WLT)
C
C        The variable NPAIRS is now initialized
C        at the same point as NA, NTHETA, ND, and NW to be
C        zero. This prevents the routine from performing
C        needless calculations for planets and avoids possible
C        floating point exceptions.
C
C-     SPICELIB Version 3.2.0, 22-MAR-1995 (KSZ)
C
C        TISBOD now gets the TSIPM matrix from PCKMAT.
C        Reference frame calculation moved to end.
C
C-     SPICELIB Version 3.0.1, 07-OCT-1994 (KSZ)
C
C        TISBOD bug which mistakenly moded DW by 2PI
C        was removed.
C
C-     SPICELIB Version 3.0.0, 10-MAR-1994 (KSZ)
C
C        TISBOD now uses new software to check for the
C        existence of binary PCK files, search the for
C        data corresponding to the requested body and time,
C        and return the appropriate Euler angles. Otherwise
C        the code calculates the Euler angles from the
C        P_constants kernel file.
C
C-     SPICELIB Version 2.0.0, 04-SEP-1991 (NJB)
C
C         Updated to handle P_constants referenced to different epochs
C         and inertial reference frames.
C
C         TISBOD now checks the kernel pool for presence of the
C         variables
C
C            BODY#_CONSTANTS_REF_FRAME
C
C         and
C
C            BODY#_CONSTANTS_JED_EPOCH
C
C         where # is the NAIF integer code of the barycenter of a
C         planetary system or of a body other than a planet or
C         satellite. If either or both of these variables are
C         present, the P_constants for BODY are presumed to be
C         referenced to the specified inertial frame or epoch.
C         If the epoch of the constants is not J2000, the input
C         time ET is converted to seconds past the reference epoch.
C         If the frame of the constants is not the frame specified
C         by REF, the rotation from the P_constants' frame to
C         body-fixed coordinates is transformed to the rotation from
C         the requested frame to body-fixed coordinates. The same
C         transformation is applied to the derivative of this
C         rotation.
C
C         Due to the prescience of the original author, the code
C         was already prepared to handle the possibility of
C         specification of a P_constants inertial reference frame via
C         kernel pool variables.
C
C
C         Also, the $Required_Reading and $Literature_References
C         sections were updated. The SPK required reading has been
C         deleted from the $Literature_References section, and the
C         NAIF_IDS, KERNEL, and TIME Required Reading files have
C         been added in the $Required_Reading section.
C
C-&

C
C     SPICELIB functions
C
      DOUBLE PRECISION      TWOPI
      DOUBLE PRECISION      HALFPI
      DOUBLE PRECISION      J2000
      DOUBLE PRECISION      RPD
      DOUBLE PRECISION      SPD
      DOUBLE PRECISION      VDOTG

      INTEGER               ZZBODBRY

      LOGICAL               BODFND
      LOGICAL               FAILED
      LOGICAL               RETURN

C
C     Local parameters
C
      INTEGER               MAXANG
      PARAMETER           ( MAXANG = 100 )

      INTEGER               MXPOLY
      PARAMETER           ( MXPOLY = 3 )

      INTEGER               NAMLEN
      PARAMETER           ( NAMLEN = 32 )

      INTEGER               FRNMLN
      PARAMETER           ( FRNMLN = 32 )

      INTEGER               TIMLEN
      PARAMETER           ( TIMLEN = 35 )

C
C     Local variables
C
      CHARACTER*(1)         DTYPE
      CHARACTER*(LMSGLN)    ERRMSG
      CHARACTER*(FRNMLN)    FIXFRM
      CHARACTER*(NAMLEN)    ITEM
      CHARACTER*(TIMLEN)    TIMSTR

      DOUBLE PRECISION      AC     (    MAXANG )
      DOUBLE PRECISION      COSTH  (    MAXANG )
      DOUBLE PRECISION      D
      DOUBLE PRECISION      DC     (    MAXANG )
      DOUBLE PRECISION      DCOEF  ( 0:MXPOLY-1 )
      DOUBLE PRECISION      DCOSTH (    MAXANG )
      DOUBLE PRECISION      DDEC
      DOUBLE PRECISION      DDELTA
      DOUBLE PRECISION      DEC
      DOUBLE PRECISION      DELTA
      DOUBLE PRECISION      DPHI
      DOUBLE PRECISION      DRA
      DOUBLE PRECISION      DSINTH (    MAXANG )
      DOUBLE PRECISION      DTIPM  (     3,  3 )
      DOUBLE PRECISION      DTHETA
      DOUBLE PRECISION      DW
      DOUBLE PRECISION      EULSTA (     6     )
      DOUBLE PRECISION      EPOCH
      DOUBLE PRECISION      PCKEPC
      DOUBLE PRECISION      PCKREF
      DOUBLE PRECISION      PHI
      DOUBLE PRECISION      RA
      DOUBLE PRECISION      RCOEF  ( 0:MXPOLY-1 )
      DOUBLE PRECISION      REQ2PC (     3, 3  )
      DOUBLE PRECISION      SINTH  (    MAXANG )
      DOUBLE PRECISION      T
      DOUBLE PRECISION      TCOEF  ( 2, MAXANG )
      DOUBLE PRECISION      TIPM   (     3, 3  )
      DOUBLE PRECISION      THETA
      DOUBLE PRECISION      TMPEPC
      DOUBLE PRECISION      W
      DOUBLE PRECISION      WC     (    MAXANG )
      DOUBLE PRECISION      WCOEF  ( 0:MXPOLY-1 )
      DOUBLE PRECISION      XDTIPM (     3, 3  )
      DOUBLE PRECISION      XTIPM  (     3, 3  )

      INTEGER               CENT
      INTEGER               DIM
      INTEGER               FRCODE
      INTEGER               I
      INTEGER               J
      INTEGER               J2CODE
      INTEGER               NA
      INTEGER               ND
      INTEGER               NPAIRS
      INTEGER               NTHETA
      INTEGER               NW
      INTEGER               PCREF
      INTEGER               REFID
      INTEGER               REQREF

      LOGICAL               FIRST
      LOGICAL               FOUND

C
C     POOL state counter.
C
      INTEGER               PULCTR ( CTRSIZ )

      LOGICAL               UPDATE

C
C     ID-based hash for text PCK data. KIDLST, KIDPOL, and KIDIDS
C     provide the index in the body PCK data arrays at which the data
C     of the body with a given ID are stored.
C
      INTEGER               LBPOOL
      PARAMETER           ( LBPOOL = -5 )

      INTEGER               MAXBOD
      PARAMETER           ( MAXBOD = 157 )

      INTEGER               BIDLST (             MAXBOD )
      INTEGER               BIDPOL ( LBPOOL :    MAXBOD )
      INTEGER               BIDIDS (             MAXBOD )

      DOUBLE PRECISION      BAC    (     MAXANG, MAXBOD )
      DOUBLE PRECISION      BDC    (     MAXANG, MAXBOD )
      DOUBLE PRECISION      BDCOEF ( 0:MXPOLY-1, MAXBOD )
      DOUBLE PRECISION      BPCKEP (             MAXBOD )
      DOUBLE PRECISION      BRCOEF ( 0:MXPOLY-1, MAXBOD )
      DOUBLE PRECISION      BTCOEF ( 2,  MAXANG, MAXBOD )
      DOUBLE PRECISION      BWC    (     MAXANG, MAXBOD )
      DOUBLE PRECISION      BWCOEF ( 0:MXPOLY-1, MAXBOD )

      INTEGER               BNA    (             MAXBOD )
      INTEGER               BND    (             MAXBOD )
      INTEGER               BNPAIR (             MAXBOD )
      INTEGER               BNW    (             MAXBOD )
      INTEGER               BPCREF (             MAXBOD )

      DOUBLE PRECISION      SINTMP
      DOUBLE PRECISION      COSTMP

      INTEGER               AT
      INTEGER               AVAIL

      LOGICAL               NEW

C
C     Saved variables
C
C     Because we need to save almost everything, we save everything
C     rather than taking a chance and accidentally leaving something
C     off the list.
C
      SAVE

C
C     Initial values
C
      DATA  FIRST  / .TRUE. /
      DATA  FOUND  / .FALSE. /

C
C     Standard SPICE error handling.
C
      IF ( RETURN () ) THEN
         RETURN
      ELSE
         CALL CHKIN ( 'TISBOD' )
      END IF

C
C     Perform any needed first pass initializations.
C
      IF ( FIRST ) THEN

C
C        Initialize POOL state counter to the user value.
C
         CALL ZZCTRUIN ( PULCTR )

C
C        Initialize kernel POOL frame hashes.
C
         CALL ZZHSIINI ( MAXBOD, BIDLST, BIDPOL )

C
C        Get the code for the J2000 frame.
C
         CALL IRFNUM ( 'J2000', J2CODE )

C
C        Set seconds per day and per century.
C
         D = SPD()
         T = D  *  36525.D0

         FIRST = .FALSE.

      END IF

      CALL IRFNUM  ( REF, REQREF )

C
C     Get state transformation matrix from high precision PCK file, if
C     available.
C
      CALL PCKMAT ( BODY, ET, PCREF, TSIPM, FOUND )

      IF ( .NOT. FOUND ) THEN

C
C        The data for the frame of interest are not available in a
C        loaded binary PCK file. This is not an error: the data may be
C        present in the kernel pool.
C
C        Check the POOL counter. If it changed, dump all buffered
C        constants data.
C
         CALL ZZPCTRCK ( PULCTR, UPDATE )

         IF ( UPDATE ) THEN
            CALL ZZHSIINI ( MAXBOD, BIDLST, BIDPOL )
         END IF

C
C        Check if we have data for this body in our buffers.
C
         CALL ZZHSICHK ( BIDLST, BIDPOL, BIDIDS, BODY, AT )

         IF ( AT .NE. 0 ) THEN

C
C           Set PCREF as it is used after the text PCK "IF".
C
            PCREF  = BPCREF ( AT )

         ELSE

C
C           Conduct a non-error-signaling check for the presence of a
C           kernel variable that is required to implement an IAU style
C           body-fixed reference frame. If the data aren't available,
C           we don't want BODVCD to signal a SPICE(KERNELVARNOTFOUND)
C           error; we want to issue the error signal locally, with a
C           better error message.
C
            ITEM = 'BODY#_PM'
            CALL REPMI  ( ITEM, '#',   BODY, ITEM  )
            CALL DTPOOL ( ITEM, FOUND, NW,   DTYPE )

            IF ( .NOT. FOUND ) THEN
C
C              Now we do have an error.
C
C              We don't have the data we'll need to produced the
C              requested state transformation matrix. In order to
C              create an error message understandable to the user,
C              find, if possible, the name of the reference frame
C              associated with the input body. Note that the body is
C              really identified by a PCK frame class ID code, though
C              most of the documentation just calls it a body ID code.
C
               CALL CCIFRM ( PCK, BODY,  FRCODE, FIXFRM, CENT, FOUND )
               CALL ETCAL  ( ET,  TIMSTR )

               ERRMSG = 'PCK data required to compute the orientation '
     .      //       'of the # # for epoch # TDB were not found. '
     .      //       'If these data were to be provided by a binary '
     .      //       'PCK file, then it is possible that the PCK '
     .      //       'file does not have coverage for the specified '
     .      //       'body-fixed frame at the time of interest. If the '
     .      //       'data were to be provided by a text PCK file, '
     .      //       'then possibly the file does not contain data '
     .      //       'for the specified body-fixed frame. In '
     .      //       'either case it is possible that a required '
     .      //       'PCK file was not loaded at all.'

C
C              Fill in the variable data in the error message.
C
               IF ( FOUND ) THEN
C
C                 The frame system knows the name of the body-fixed
C                 frame.
C
                  CALL SETMSG ( ERRMSG                  )
                  CALL ERRCH  ( '#', 'body-fixed frame' )
                  CALL ERRCH  ( '#', FIXFRM             )
                  CALL ERRCH  ( '#', TIMSTR             )

               ELSE
C
C                 The frame system doesn't know the name of the
C                 body-fixed frame, most likely due to a missing
C                 frame kernel.
C
                  CALL SUFFIX ( '#', 1, ERRMSG                      )
                  CALL SETMSG ( ERRMSG                              )
                  CALL ERRCH  ( '#',
     .                       'body-fixed frame associated with '
     .         //            'the ID code'                       )
                  CALL ERRINT ( '#', BODY                           )
                  CALL ERRCH  ( '#', TIMSTR                         )
                  CALL ERRCH  ( '#',
     .                       'Also, a frame kernel defining the '
     .         //            'body-fixed frame associated with '
     .         //            'body # may need to be loaded.'     )
                  CALL ERRINT ( '#', BODY                           )

               END IF

               CALL SIGERR ( 'SPICE(FRAMEDATANOTFOUND)' )
               CALL CHKOUT ( 'TISBOD'                   )
               RETURN

            END IF

C
C           Find the body code used to label the reference frame and
C           epoch specifiers for the orientation constants for BODY.
C
C           For planetary systems, the reference frame and epoch for
C           the orientation constants is associated with the system
C           barycenter, not with individual bodies in the system. For
C           any other bodies, (the Sun or asteroids, for example) the
C           body's own code is used as the label.
C
            REFID = ZZBODBRY ( BODY )

C
C           Look up the epoch of the constants. The epoch is specified
C           as a Julian ephemeris date. The epoch defaults to J2000.
C
            ITEM  =  'BODY#_CONSTANTS_JED_EPOCH'

            CALL REPMI  ( ITEM, '#', REFID, ITEM               )
            CALL GDPOOL ( ITEM,  1,  1,     DIM, PCKEPC, FOUND )

            IF ( .NOT. FOUND ) THEN

               PCKEPC = J2000()

            END IF

C
C           Look up the reference frame of the constants. The reference
C           frame is specified by a code recognized by CHGIRF. The
C           default frame is J2000, symbolized by the code J2CODE.
C
            ITEM  =  'BODY#_CONSTANTS_REF_FRAME'

            CALL REPMI  ( ITEM, '#', REFID, ITEM               )
            CALL GDPOOL ( ITEM,  1,  1,     DIM, PCKREF, FOUND )

            IF ( FOUND ) THEN
               PCREF  =  NINT ( PCKREF )
            ELSE
               PCREF  =  J2CODE
            END IF

C
C           Whatever the body, it has quadratic time polynomials for
C           the RA and Dec of the pole, and for the rotation of the
C           Prime Meridian.
C
            ITEM = 'POLE_RA'
            CALL CLEARD (             MXPOLY,     RCOEF )
            CALL BODVCD ( BODY, ITEM, MXPOLY, NA, RCOEF )

            IF ( FAILED() ) THEN
               CALL CHKOUT ( 'TISBOD' )
               RETURN
            END IF

            ITEM = 'POLE_DEC'
            CALL CLEARD (             MXPOLY,     DCOEF )
            CALL BODVCD ( BODY, ITEM, MXPOLY, ND, DCOEF )

            IF ( FAILED() ) THEN
               CALL CHKOUT ( 'TISBOD' )
               RETURN
            END IF

            ITEM = 'PM'
            CALL CLEARD (             MXPOLY,     WCOEF )
            CALL BODVCD ( BODY, ITEM, MXPOLY, NW, WCOEF )

            IF ( FAILED() ) THEN
               CALL CHKOUT ( 'TISBOD' )
               RETURN
            END IF

C
C           If the body is a satellite, there may be additional
C           nutation and libration (THETA) terms.
C
            NTHETA = 0
            NPAIRS = 0
            NA     = 0
            ND     = 0
            NW     = 0

            ITEM = 'NUT_PREC_ANGLES'

C
C           There is something a bit obscure going on below. We are
C           passing a two dimensional array ( TCOEF(2, MAXANG) ). But
C           BODVCD is expecting a 1- dimensional array. BODVCD loads
C           the array TCOEF in the following order
C
C              TCOEF(1,1), TCOEF(2,1), TCOEF(1,2), TCOEF(2,2),   ...
C
C           The NTHETA that comes back is the total number of items
C           loaded, but we will need the actual limit on the second
C           dimension. That is --- NTHETA / 2.
C
            IF ( BODFND ( REFID, ITEM ) ) THEN

               CALL BODVCD ( REFID, ITEM, MAXANG*2, NTHETA, TCOEF )

               IF ( FAILED() ) THEN
                  CALL CHKOUT ( 'TISBOD' )
                  RETURN
               END IF

               NPAIRS = NTHETA / 2

            END IF

C
C           Look up the right ascension nutations in the precession of
C           the pole. NA is the number of Ascension coefficients. AC
C           are the Ascension coefficients.
C
            ITEM = 'NUT_PREC_RA'

            IF ( BODFND ( BODY, ITEM ) ) THEN

               CALL BODVCD ( BODY, ITEM, MAXANG, NA, AC )

               IF ( FAILED() ) THEN
                  CALL CHKOUT ( 'TISBOD' )
                  RETURN
               END IF

            END IF

C
C           Look up the declination nutations in the precession of the
C           pole. ND is the number of Declination coefficients. DC are
C           the Declination coefficients.
C
            ITEM = 'NUT_PREC_DEC'

            IF ( BODFND ( BODY, ITEM ) ) THEN

               CALL BODVCD ( BODY, ITEM, MAXANG, ND, DC )

               IF ( FAILED() ) THEN
                  CALL CHKOUT ( 'TISBOD' )
                  RETURN
               END IF

            END IF

C
C           Finally look up the prime meridian nutations. NW is the
C           number of coefficients. WC is the array of coefficients.
C
            ITEM = 'NUT_PREC_PM'

            IF ( BODFND ( BODY, ITEM ) ) THEN

               CALL BODVCD ( BODY, ITEM, MAXANG, NW, WC )

               IF ( FAILED() ) THEN
                  CALL CHKOUT ( 'TISBOD' )
                  RETURN
               END IF

            END IF

C
C           The number of coefficients returned had better not be
C           bigger than the number of angles we are going to compute.
C           If it is we simply signal an error and bag it, for sure.
C
            IF ( MAX ( NA, ND, NW ) .GT. NPAIRS ) THEN
               CALL SETMSG ( 'Insufficient number of nutation/' //
     .                    'precession angles for body * at time #.' )
               CALL ERRINT ( '*', BODY )
               CALL ERRDP  ( '#', ET   )
               CALL SIGERR ( 'SPICE(INSUFFICIENTANGLES)' )
               CALL CHKOUT ( 'TISBOD' )
               RETURN
            END IF

C
C           We succeeded to fetch all data for this body. Save it in
C           the buffers. First check if there is room. If not, dump the
C           buffers to make room.
C
            CALL ZZHSIAVL( BIDPOL, AVAIL )
            IF ( AVAIL .LE. 0 ) THEN
               CALL ZZHSIINI ( MAXBOD, BIDLST, BIDPOL )
            END IF

C
C           Add this body to the hash and save its data in the buffers.
C
            CALL ZZHSIADD ( BIDLST, BIDPOL, BIDIDS, BODY, AT, NEW )

            BPCKEP ( AT ) = PCKEPC
            BPCREF ( AT ) = PCREF
            BNPAIR ( AT ) = NPAIRS
            BNA    ( AT ) = NA
            BND    ( AT ) = ND
            BNW    ( AT ) = NW

            CALL MOVED ( RCOEF ( 0 ),    MXPOLY,   BRCOEF ( 0,    AT ) )
            CALL MOVED ( DCOEF ( 0 ),    MXPOLY,   BDCOEF ( 0,    AT ) )
            CALL MOVED ( WCOEF ( 0 ),    MXPOLY,   BWCOEF ( 0,    AT ) )
            CALL MOVED ( TCOEF ( 1, 1 ), MAXANG*2, BTCOEF ( 1, 1, AT ) )
            CALL MOVED ( AC    ( 1 ),    MAXANG,   BAC    ( 1,    AT ) )
            CALL MOVED ( DC    ( 1 ),    MAXANG,   BDC    ( 1,    AT ) )
            CALL MOVED ( WC    ( 1 ),    MAXANG,   BWC    ( 1,    AT ) )

         END IF

C
C        The reference epoch in the PCK is given as JED. Convert to
C        ephemeris seconds past J2000. Then convert the input ET to
C        seconds past the reference epoch.
C
         TMPEPC  =  SPD() * ( BPCKEP(AT) - J2000() )
         EPOCH   =  ET    -   TMPEPC

C
C        Evaluate the time polynomials and their derivatives w.r.t.
C        EPOCH at EPOCH.
C
C        Evaluate the time polynomials at EPOCH.
C
         RA  = BRCOEF(0,AT) +
     .         (EPOCH/T) * ( BRCOEF(1,AT) + (EPOCH/T) * BRCOEF(2,AT) )
         DEC = BDCOEF(0,AT) +
     .         (EPOCH/T) * ( BDCOEF(1,AT) + (EPOCH/T) * BDCOEF(2,AT) )
         W   = BWCOEF(0,AT) +
     .         (EPOCH/D) * ( BWCOEF(1,AT) + (EPOCH/D) * BWCOEF(2,AT) )

         DRA  = ( BRCOEF(1,AT)  + 2.0D0*(EPOCH/T)*BRCOEF(2,AT) ) / T
         DDEC = ( BDCOEF(1,AT)  + 2.0D0*(EPOCH/T)*BDCOEF(2,AT) ) / T
         DW   = ( BWCOEF(1,AT)  + 2.0D0*(EPOCH/D)*BWCOEF(2,AT) ) / D

C
C        Compute the nutations and librations (and their derivatives)
C        as appropriate.
C
         DO   I = 1, BNPAIR(AT)

            THETA  = ( BTCOEF(1,I,AT) + (EPOCH/T) * BTCOEF(2,I,AT) )
     .                                                         * RPD()
            DTHETA = ( BTCOEF(2,I,AT) / T ) * RPD()

            SINTMP    =  SIN ( THETA )
            COSTMP    =  COS ( THETA )

            SINTH(I)  =  SINTMP
            COSTH(I)  =  COSTMP
            DSINTH(I) =  COSTMP * DTHETA
            DCOSTH(I) = -SINTMP * DTHETA

         END DO

C
C        Adjust RA, DEC, W and their derivatives by the librations
C        and nutations.
C
         RA  = RA    + VDOTG ( BAC(1,AT), SINTH,  BNA(AT) )
         DEC = DEC   + VDOTG ( BDC(1,AT), COSTH,  BND(AT) )
         W   = W     + VDOTG ( BWC(1,AT), SINTH,  BNW(AT) )

         DRA  = DRA  + VDOTG ( BAC(1,AT), DSINTH, BNA(AT) )
         DDEC = DDEC + VDOTG ( BDC(1,AT), DCOSTH, BND(AT) )
         DW   = DW   + VDOTG ( BWC(1,AT), DSINTH, BNW(AT) )

C
C        Convert from degrees to radians
C
         RA      = RA   * RPD()
         DEC     = DEC  * RPD()
         W       = W    * RPD()

         DRA     = DRA  * RPD()
         DDEC    = DDEC * RPD()
         DW      = DW   * RPD()
C
C        Convert to Euler angles.
C
         W       = MOD ( W,  TWOPI() )

         PHI     = RA + HALFPI()
         DELTA   = HALFPI() - DEC

         DPHI    =   DRA
         DDELTA  = - DDEC

         IF ( FAILED() ) THEN
            CALL CHKOUT ( 'TISBOD' )
            RETURN
         END IF

C
C        Pack the Euler angles and their derivatives into
C        a state vector.
C
         CALL VPACK ( W,  DELTA,  PHI,  EULSTA    )
         CALL VPACK ( DW, DDELTA, DPHI, EULSTA(4) )

C
C        Find the state transformation defined by the Euler angle
C        state vector. The transformation matrix TSIPM has the
C        following structure:
C
C            -            -
C           |       :      |
C           | TIPM  :  0   |
C           | ......:......|
C           |       :      |
C           | DTIPM : TIPM |
C           |       :      |
C            -            -
C
         CALL EUL2XF ( EULSTA, 3, 1, 3, TSIPM )

      END IF

C
C     At this point the base frame PCREF has been determined.
C
C     If the requested base frame is not base frame associated with the
C     PCK data, adjust the transformation matrix TSIPM to map from the
C     requested frame to the body-fixed frame.
C
      IF ( REQREF .NE. PCREF ) THEN
C
C        Next get the position transformation from the user specified
C        inertial frame to the native PCK inertial frame.
C
         CALL IRFROT ( REQREF, PCREF, REQ2PC )

         IF ( FAILED() ) THEN
            CALL CHKOUT ( 'TISBOD' )
            RETURN
         END IF

C
C        Since we're applying an inertial transformation to TSIPM,
C        we can rotate the non-zero blocks of TSIPM. This saves
C        a bunch of double precision multiplications.
C
C        Extract the upper and lower left blocks of TSIPM.
C
         DO I = 1, 3

            DO J = 1, 3

               TIPM (I,J) =  TSIPM( I,   J )
               DTIPM(I,J) =  TSIPM( I+3, J )

            END DO

         END DO

C
C        Rotate the blocks. Note this is a right multiplication.
C
         CALL MXM ( TIPM,  REQ2PC, XTIPM  )
         CALL MXM ( DTIPM, REQ2PC, XDTIPM )

C
C        Replace the non-zero blocks of TSIPM. This gives us the
C        transformation from the requested frame to the
C        bodyfixed frame.
C
         DO I = 1, 3

            DO J = 1, 3

               TSIPM( I,   J   ) = XTIPM (I,J)
               TSIPM( I+3, J+3 ) = XTIPM (I,J)
               TSIPM( I+3, J   ) = XDTIPM(I,J)

            END DO

         END DO

      END IF

C
C     That's all folks. Check out and get out.
C

      CALL CHKOUT ( 'TISBOD' )
      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


text:
      COMMON/ModFiles/filein,l_file,fileout1,fileout2,fileout3,
     1 l_par,filepar,l_lpt,lptfile,opt,filenorm
      integer l_file,l_par,l_lpt,opt
      character*140 filein,fileout1,fileout2,fileout3,filepar,
     1 lptfile,filenorm

max_stars_repo_name: mantidproject/quasielasticbayes

max_issues_repo_name: mantidproject/quasielasticbayes

max_forks_repo_name: mantidproject/quasielasticbayes


text:
! Copyright (c) 2019, NVIDIA CORPORATION.  All rights reserved.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.
!

module mod
logical expect(1), rslt(1)
type obj
contains
procedure :: bar
end type

contains
subroutine foo()
!print *, 'in foo'
end subroutine

subroutine foo2(i)
integer :: i
!print *, 'in foo2', i
rslt(1) = i .eq. -99
end subroutine

subroutine bar(this)
class(obj) :: this
block
type t
contains
procedure, nopass :: foo
procedure, nopass :: foo2
generic :: func => foo, foo2
end type
type(t) :: o

call o%func(-99)
end block
end subroutine

end module

use mod
type(obj) :: x
expect = .true.
rslt = .false.
call x%bar()
call check(rslt, expect, 1)
end

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
submodule (meshing) meshing_operators

  implicit none

contains

  !>  Returns the distance between two cell centres
  module subroutine get_neighbour_distance(loc_p, loc_nb, dx)
    type(cell_locator), intent(in) :: loc_p             !< The cell distance is measured from.
    type(neighbour_locator), intent(in) :: loc_nb       !< The cell distance is measured to.
    real(ccs_real), dimension(ndim), intent(out) :: dx  !< ndim-array of the distance

    real(ccs_real), dimension(ndim) :: xp
    real(ccs_real), dimension(ndim) :: xnb

    call get_centre(loc_p, xp)
    call get_centre(loc_nb, xnb)

    dx(:) = xnb(:) - xp(:)

  end subroutine get_neighbour_distance

  !>  Returns the distance from cell to face centres
  module subroutine get_face_distance(loc_p, loc_f, dx)
    type(cell_locator), intent(in) :: loc_p             !< The cell distance is measured from.
    type(face_locator), intent(in) :: loc_f             !< The face distance is measured to.
    real(ccs_real), dimension(ndim), intent(out) :: dx  !< ndim-array of the distance

    real(ccs_real), dimension(ndim) :: xp
    real(ccs_real), dimension(ndim) :: xf

    call get_centre(loc_p, xp)
    call get_centre(loc_f, xf)

    dx(:) = xf(:) - xp(:)

  end subroutine get_face_distance

end submodule meshing_operators

max_stars_repo_name: asimovpp/asimov-ccs

max_issues_repo_name: asimovpp/asimov-ccs

max_forks_repo_name: asimovpp/asimov-ccs


! { dg-do compile }
! Tests the fix for a further regression caused by the
! fix for PR28788 and posted as PR28908. The problem was
! caused by the patch preventing interface derived types
! from associating with identical derived types in the
! containing namespaces.
!
! Contributed by HJ Lu  <hjl@lucon.org>
!
module bar
  implicit none
  public
  type domain_ptr
    type(domain), POINTER  :: ptr
  end type domain_ptr
  type domain
    TYPE(domain_ptr) , DIMENSION( : ) , POINTER         :: parents
    TYPE(domain_ptr) , DIMENSION( : ) , POINTER         :: nests
  end type domain
end module bar

module foo
contains
  recursive subroutine integrate (grid)
    use bar
    implicit none
    type(domain), POINTER  :: grid
    interface
      subroutine solve_interface (grid)
        use bar
        TYPE (domain) grid
      end subroutine solve_interface
    end interface
  end subroutine integrate
end module foo


!*****************************************************************************************
!> author: Jacob Williams
!
!  Define the numeric kinds.

    module kind_module

    use, intrinsic :: iso_fortran_env,    only: real32,real64,real128

    implicit none

    private

    !integer,parameter,public :: wp = real32     !! single precision reals
    integer,parameter,public :: wp = real64      !! double precision reals
    !integer,parameter,public :: wp = real128    !! quad precision reals

    end module kind_module
!*****************************************************************************************