text:
c Author: Anand Patil, anand.prabhakar.patil@gmail.com

      SUBROUTINE checksymm(X,n,cs)
c Checks symmetry
cf2py intent(hide) n
cf2py intent(out) cs
      DOUBLE PRECISION X(n,n)
      INTEGER n,i,j
      LOGICAL cs

      cs = .FALSE.
      do j=1,n-1
          do i=j+1,n
              if (X(i,j).NE.X(j,i)) then
                  cs = .TRUE.
                  return
              end if
          end do
      end do

      return
      END

      SUBROUTINE chol_mvnorm(x, mu, sig, n, like, info)

cf2py double precision dimension(n), intent(copy) :: x
cf2py double precision dimension(n), intent(copy) :: mu
cf2py integer intent(hide),depend(x) :: n=len(x)
cf2py double precision dimension(n,n), intent(in) :: sig
cf2py double precision intent(out) :: like
cf2py integer intent(hide) :: info
cf2py threadsafe

      DOUBLE PRECISION sig(n,n), x(n), mu(n), like
      INTEGER n, info, i
      DOUBLE PRECISION twopi_N, log_detC, gd
      DOUBLE PRECISION infinity
      PARAMETER (infinity = 1.7976931348623157d308)
      DOUBLE PRECISION PI
      PARAMETER (PI=3.141592653589793238462643d0)

      EXTERNAL DPOTRS
! DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO ) Solves triangular system
      EXTERNAL DAXPY
! DAXPY(N,DA,DX,INCX,DY,INCY) Adding vectors
      EXTERNAL DCOPY
! DCOPY(N,DX,INCX,DY,INCY) copies x to y
! NB DDOT from ATLAS, compiled with gfortran 4.2 on Ubuntu Gutsy,
! was producing bad output- hence the manual dot product.

!     x <- (x-mu)
      call DAXPY(n, -1.0D0, mu, 1, x, 1)

!       mu <- x
      call DCOPY(n,x,1,mu,1)

!     x <- sig ^-1 * x
      call DPOTRS('L',n,1,sig,n,x,n,info)

      gd=0.0D0
      do i=1,n
          gd=gd+x(i)*mu(i)
      end do

!     like <- .5 dot(x,mu) (.5 (x-mu) C^{-1} (x-mu)^T)
      like = -0.5D0 * gd
!       print *, like

      twopi_N = 0.5D0 * N * dlog(2.0D0*PI)
!       print *, twopi_N

      log_detC = 0.0D0
      do i=1,n
        log_detC = log_detC + log(sig(i,i))
      enddo
!       print *, log_detC

      like = like - twopi_N - log_detC

      return
      END

      SUBROUTINE cov_mvnorm(x, mu, C, n, like, info)

cf2py double precision dimension(n), intent(copy) :: x
cf2py double precision dimension(n), intent(copy) :: mu
cf2py integer intent(hide),depend(x) :: n=len(x)
cf2py double precision dimension(n,n), intent(copy) :: C
cf2py double precision intent(out) :: like
cf2py integer intent(hide) :: info
cf2py threadsafe

      DOUBLE PRECISION C(n,n), x(n), mu(n), like
      INTEGER n, info
      DOUBLE PRECISION infinity
      PARAMETER (infinity = 1.7976931348623157d308)
      DOUBLE PRECISION PI
      PARAMETER (PI=3.141592653589793238462643d0)

      EXTERNAL DPOTRF
! DPOTRF( UPLO, N, A, LDA, INFO ) Cholesky factorization

!     C <- cholesky(C)
      call DPOTRF( 'L', n, C, n, info )
!       print *, C

!      If cholesky failed, puke.
       if (info .GT. 0) then
         like = -infinity
         RETURN
       endif

!     Call to chol_mvnorm
      call chol_mvnorm(x,mu,C,n,like,info)

      return
      END

      SUBROUTINE prec_mvnorm(x, mu, tau, n, like)

cf2py double precision dimension(n), intent(copy) :: x
cf2py double precision dimension(n), intent(copy) :: mu
cf2py integer intent(hide),depend(x) :: n=len(x)
cf2py double precision dimension(n,n), intent(copy) :: tau
cf2py double precision intent(out) :: like
cf2py threadsafe

      DOUBLE PRECISION tau(n,n), x(n), mu(n), like, gd
      INTEGER n, info, i
      DOUBLE PRECISION infinity
      PARAMETER (infinity = 1.7976931348623157d308)
      DOUBLE PRECISION PI
      PARAMETER (PI=3.141592653589793238462643d0)
      DOUBLE PRECISION twopi_N, log_dettau

      EXTERNAL DPOTRF
! DPOTRF( UPLO, N, A, LDA, INFO ) Cholesky factorization
      EXTERNAL DSYMV
! DSYMV(UPLO,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY) Symmetric matrix-vector multiply
      EXTERNAL DAXPY
! DAXPY(N,DA,DX,INCX,DY,INCY) Adding vectors
      EXTERNAL DCOPY
! DCOPY(N,DX,INCX,DY,INCY) copies x to y

      twopi_N = 0.5D0 * N * dlog(2.0D0*PI)

!     x <- (x-mu)
      call DAXPY(n, -1.0D0, mu, 1, x, 1)

!       mu <- x
      call DCOPY(n,x,1,mu,1)

!     mu <- tau * x
      call DSYMV('L',n,1.0D0,tau,n,x,1,0.0D0,mu,1)

      gd = 0.0D0
      do i=1,n
          gd=gd+x(i)*mu(i)
      end do

!     like <- -.5 dot(x,mu) (.5 (x-mu) C^{-1} (x-mu)^T)
      like = -0.5D0 * gd

!      Cholesky factorize tau for the determinant.
       call DPOTRF( 'L', n, tau, n, info )

!      If cholesky failed, puke.
       if (info .GT. 0) then
         like = -infinity
         RETURN
       endif

!      Otherwise read off determinant.
       log_dettau = 0.0D0
       do i=1,n
         log_dettau = log_dettau + dlog(tau(i,i))
       enddo

       like = like - twopi_N + log_dettau

      return
      END

      SUBROUTINE blas_wishart(X,k,n,T,like)

c Wishart log-likelihood function.

cf2py double precision dimension(k,k),intent(copy) :: X,T
cf2py double precision intent(out) :: like
cf2py integer intent(hide),depend(X) :: k=len(X)
cf2py threadsafe

      INTEGER i,k,n
      DOUBLE PRECISION X(k,k),T(k,k),bx(k,k)
      DOUBLE PRECISION dx,db,tbx,a,g,like
      DOUBLE PRECISION infinity
      PARAMETER (infinity = 1.7976931348623157d308)
      DOUBLE PRECISION PI
      LOGICAL csx
      PARAMETER (PI=3.141592653589793238462643d0)

      EXTERNAL DCOPY
! DCOPY(N,DX,INCX,DY,INCY) copies x to y
      EXTERNAL DSYMM
! DSYMM(SIDE,UPLO,M,N,ALPHA,A,LDA,B,LDB,BETA,C,LDC) alpha*A*B + beta*C when side='l'
      EXTERNAL DPOTRF
! DPOTRF( UPLO, N, A, LDA, INFO ) Cholesky factorization

c Check X for symmetry
      CALL checksymm(X,k,csx)
      if (csx) then
          like = -infinity
          return
      end if

c trace of T*X
!     bx <- T * X
      call DSYMM('L','L',k,k,1.0D0,T,k,x,k,0.0D0,bx,k)

c Cholesky factor T, puke if not pos def.
      call DPOTRF( 'L', k, T, k, info )
      if (info .GT. 0) then
        like = -infinity
        RETURN
      endif

c Cholesky factor X, puke if not pos def.
      call DPOTRF( 'L', k, X, k, info )
      if (info .GT. 0) then
        like = -infinity
        RETURN
      endif

c Get the trace and log-sqrt-determinants
      tbx = 0.0D0
      dx = 0.0D0
      db = 0.0D0

      do i=1,k
        tbx = tbx + bx(i,i)
        dx = dx + dlog(X(i,i))
        db = db + dlog(T(i,i))
      enddo

      if (k .GT. n) then
        like = -infinity
        RETURN
      endif

      like = 0.0D0
      like = (n - k - 1) * dx
      like = like + n * db
      like = like - 0.5D0 * tbx
      like = like - (n*k/2.0d0)*dlog(2.0d0)

      do i=1,k
        a = (n - i + 1)/2.0D0
        call gamfun(a, g)
        like = like - g
      enddo

      like = like - k * (k-1) * 0.25D0 * dlog(PI)
!
      return
      END

c

      SUBROUTINE blas_wishart_cov(X,k,n,V,like)

c Wishart log-likelihood function.
c Doesn't vectorize the determinants, just the matrix multiplication.

cf2py double precision dimension(k,k),intent(copy) :: V,X
cf2py double precision intent(out) :: like
cf2py integer intent(hide),depend(X) :: k=len(X)
cf2py threadsafe

      INTEGER i,k,info, n
      DOUBLE PRECISION X(k,k),V(k,k),bx(k,k)
      DOUBLE PRECISION dx,db,tbx,a,g,like
      DOUBLE PRECISION infinity
      PARAMETER (infinity = 1.7976931348623157d308)
      DOUBLE PRECISION PI
      LOGICAL csx
      PARAMETER (PI=3.141592653589793238462643d0)
c
      EXTERNAL DCOPY
! DCOPY(N,DX,INCX,DY,INCY) copies x to y
      EXTERNAL DPOTRF
! DPOTRF( UPLO, N, A, LDA, INFO ) Cholesky factorization
      EXTERNAL DPOTRS
! DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO ) Solves triangular system

c determinants
c Check X for symmetry
      CALL checksymm(X,k,csx)
      if (csx) then
          like = -infinity
          return
      end if

c Cholesky factorize sigma, puke if not pos def
!     V <- cholesky(V)
      call DPOTRF( 'L', k, V, k, info )
      if (info .GT. 0) then
        like = -infinity
        RETURN
      endif

c trace of V^{-1}*X
!     bx <- X
      call DCOPY(k * k,X,1,bx,1)
!     bx <- V^{-1} * bx
      call DPOTRS('L',k,k,V,k,bx,k,info)

!     X <- cholesky(X)
      call DPOTRF( 'L', k, X, k, info )

! sqrt-log-determinant of sigma and X, and trace
      db=0.0D0
      dx=0.0D0
      tbx = 0.0D0
      do i=1,k
        db = db + dlog(V(i,i))
        dx = dx + dlog(X(i,i))
        tbx = tbx + bx(i,i)
      enddo

      if (k .GT. n) then
        like = -infinity
        RETURN
      endif

      like = (n - k - 1) * dx
      like = like - n * db
      like = like - 0.5D0*tbx
      like = like - (n*k/2.0D0)*dlog(2.0d0)

      do i=1,k
        a = (n - i + 1)/2.0D0
        call gamfun(a, g)
        like = like - g
      enddo

      like = like - k * (k-1) * 0.25D0 * dlog(PI)

      return
      END

      SUBROUTINE blas_inv_wishart(X,k,n,T,like)

c Inverse Wishart log-likelihood function.

cf2py double precision dimension(k,k),intent(copy) :: X,T
cf2py double precision intent(out) :: like
cf2py integer intent(hide),depend(X) :: k=len(X)
cf2py threadsafe

      INTEGER i,k,n
      DOUBLE PRECISION X(k,k),T(k,k),TX(k,k)
      DOUBLE PRECISION dx,db,tbx,a,g,like
      DOUBLE PRECISION infinity
      PARAMETER (infinity = 1.7976931348623157d308)
      DOUBLE PRECISION PI
      LOGICAL csx
      PARAMETER (PI=3.141592653589793238462643d0)

      EXTERNAL DCOPY
! DCOPY(N,DX,INCX,DY,INCY) copies x to y
      EXTERNAL DPOTRF
! DPOTRF( UPLO, N, A, LDA, INFO ) Cholesky factorization
      EXTERNAL DTRMM
! DTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB) B := alpha*B*op( A )

c Check X and T for symmetry
      CALL checksymm(X,k,csx)
      if (csx) then
          like = -infinity
          return
      end if

c trace of T*X^{-1}
!     TX <- T
      call DCOPY(k*k,T,1,TX,1)
!     TX <- T * X^{-1}
      call DTRMM('R','L','T','N',k,k,1.0D0,X,k,TX,k)

c Cholesky factor T, puke if not pos def.
      call DPOTRF( 'L', k, T, k, info )
      if (info .GT. 0) then
        like = -infinity
        RETURN
      endif

c Cholesky factor X, puke if not pos def.
      call DPOTRF( 'L', k, X, k, info )
      if (info .GT. 0) then
        like = -infinity
        RETURN
      endif

c Get the trace and log-sqrt-determinants
      tbx = 0.0D0
      dx = 0.0D0
      db = 0.0D0

      do i=1,k
        tbx = tbx + TX(i,i)
        dx = dx + dlog(X(i,i))
        db = db + dlog(T(i,i))
      enddo

      if (k .GT. n) then
        like = -infinity
        RETURN
      endif

      like = -1.0D0*(n + k + 1) * dx
      like = like + n * db
      like = like - 0.5D0 * tbx
      like = like - (n*k/2.0d0)*dlog(2.0d0)

      do i=1,k
        a = (n - i + 1)/2.0D0
        call gamfun(a, g)
        like = like - g
      enddo

      like = like - k * (k-1) * 0.25D0 * dlog(PI)
!
      return
      END

      SUBROUTINE dtrsm_wrap(M,N,A,B,SIDE,TRANSA,UPLO)
cf2py intent(inplace)::B
cf2py integer intent(hide)::M
cf2py integer intent(hide)::N
cf2py threadsafe

*     .. Scalar Arguments ..
      DOUBLE PRECISION ALPHA
      INTEGER LDA,LDB,M,N
      CHARACTER DIAG,SIDE,TRANSA,UPLO
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION A(M,M),B(M,N)

      EXTERNAL DTRSM
! DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
*     op( A )*X = alpha*B,   or   X*op( A ) = alpha*B,

      DIAG = 'N'
      LDA = M
      LDB = M
      ALPHA=1.0D0

      CALL DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
      RETURN
      END

      SUBROUTINE dtrmm_wrap(M,N,A,B,SIDE,TRANSA,UPLO)
cf2py intent(inplace)::B
cf2py integer intent(hide)::M
cf2py integer intent(hide)::N
cf2py threadsafe

*     .. Scalar Arguments ..
      DOUBLE PRECISION ALPHA
      INTEGER LDA,LDB,M,N
      CHARACTER DIAG,SIDE,TRANSA,UPLO
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION A(M,M),B(M,N)

      EXTERNAL DTRMM
! DTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
*     X = alpha*op( A )*B,   or   X = alpha*B*op( A ),

      DIAG = 'N'
      LDA = M
      LDB = M
      ALPHA=1.0D0

      CALL DTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
      RETURN
      END

      subroutine dchdc_wrap(a,p,work,piv,info)

cf2py intent(inplace) a
cf2py intent(hide) p
cf2py intent(hide) work
cf2py intent(out) piv
cf2py intent(out) info
cf2py threadsafe

      integer lda,p,piv(p),job,info
      double precision a(p,p),work(p)

c
c     dchdc computes the cholesky decomposition of a positive definite
c     matrix.  a pivoting option allows the user to estimate the
c     condition of a positive definite matrix or determine the rank
c     of a positive semidefinite matrix.
c
c     on entry
c
c         a      double precision(lda,p).
c                a contains the matrix whose decomposition is to
c                be computed.  onlt the upper half of a need be stored.
c                the lower part of the array a is not referenced.
c
c         lda    integer.
c                lda is the leading dimension of the array a.
c
c         Since f2py will make sure that a is column-major contiguous,
c         lda is set to p in this wrapper.
c
c         p      integer.
c                p is the order of the matrix.
c
c         work   double precision.
c                work is a work array.
c
c         piv   integer(p).
c                piv contains integers that control the selection
c                of the pivot elements, if pivoting has been requested.
c                each diagonal element a(k,k)
c                is placed in one of three classes according to the
c                value of piv(k).
c
c                   if piv(k) .gt. 0, then x(k) is an initial
c                                      element.
c
c                   if piv(k) .eq. 0, then x(k) is a free element.
c
c                   if piv(k) .lt. 0, then x(k) is a final element.
c
c                before the decomposition is computed, initial elements
c                are moved by symmetric row and column interchanges to
c                the beginning of the array a and final
c                elements to the end.  both initial and final elements
c                are frozen in place during the computation and only
c                free elements are moved.  at the k-th stage of the
c                reduction, if a(k,k) is occupied by a free element
c                it is interchanged with the largest free element
c                a(l,l) with l .ge. k.  piv is not referenced if
c                job .eq. 0.
c
c        job     integer.
c                job is an integer that initiates column pivoting.
c                if job .eq. 0, no pivoting is done.
c                if job .ne. 0, pivoting is done.
c
c        In this wrapper, job=1.
c
c     on return
c
c         a      a contains in its upper half the cholesky factor
c                of the matrix a as it has been permuted by pivoting.
c
c         piv   piv(j) contains the index of the diagonal element
c                of a that was moved into the j-th position,
c                provided pivoting was requested.
c
c         info   contains the index of the last positive diagonal
c                element of the cholesky factor.
c
c     for positive definite matrices info = p is the normal return.
c     for pivoting with positive semidefinite matrices info will
c     in general be less than p.  however, info may be greater than
c     the rank of a, since rounding error can cause an otherwise zero
c     element to be positive. indefinite systems will always cause
c     info to be less than p.
c
c     linpack. this version dated 08/14/78 .
c     j.j. dongarra and g.w. stewart, argonne national laboratory and
c     university of maryland.
c
c
c     blas daxpy,dswap
c     fortran dsqrt
c
c     internal variables
c
      integer pu,pl,plp1,i,j,jp,jt,k,kb,km1,kp1,l,maxl
      double precision temp
      double precision maxdia
      logical swapk,negk
c
      pl = 1
      pu = 0
      info = p
      job=1
      lda=p
      if (job .eq. 0) go to 160
c
c        pivoting has been requested. rearrange the
c        the elements according to piv.
c
         do 70 k = 1, p
            swapk = piv(k) .gt. 0
            negk = piv(k) .lt. 0
            piv(k) = k
            if (negk) piv(k) = -piv(k)
            if (.not.swapk) go to 60
               if (k .eq. pl) go to 50
                  call dswap(pl-1,a(1,k),1,a(1,pl),1)
                  temp = a(k,k)
                  a(k,k) = a(pl,pl)
                  a(pl,pl) = temp
                  plp1 = pl + 1
                  if (p .lt. plp1) go to 40
                  do 30 j = plp1, p
                     if (j .ge. k) go to 10
                        temp = a(pl,j)
                        a(pl,j) = a(j,k)
                        a(j,k) = temp
                     go to 20
   10                continue
                     if (j .eq. k) go to 20
                        temp = a(k,j)
                        a(k,j) = a(pl,j)
                        a(pl,j) = temp
   20                continue
   30             continue
   40             continue
                  piv(k) = piv(pl)
                  piv(pl) = k
   50          continue
               pl = pl + 1
   60       continue
   70    continue
         pu = p
         if (p .lt. pl) go to 150
         do 140 kb = pl, p
            k = p - kb + pl
            if (piv(k) .ge. 0) go to 130
               piv(k) = -piv(k)
               if (pu .eq. k) go to 120
                  call dswap(k-1,a(1,k),1,a(1,pu),1)
                  temp = a(k,k)
                  a(k,k) = a(pu,pu)
                  a(pu,pu) = temp
                  kp1 = k + 1
                  if (p .lt. kp1) go to 110
                  do 100 j = kp1, p
                     if (j .ge. pu) go to 80
                        temp = a(k,j)
                        a(k,j) = a(j,pu)
                        a(j,pu) = temp
                     go to 90
   80                continue
                     if (j .eq. pu) go to 90
                        temp = a(k,j)
                        a(k,j) = a(pu,j)
                        a(pu,j) = temp
   90                continue
  100             continue
  110             continue
                  jt = piv(k)
                  piv(k) = piv(pu)
                  piv(pu) = jt
  120          continue
               pu = pu - 1
  130       continue
  140    continue
  150    continue
  160 continue
      do 270 k = 1, p
c
c        reduction loop.
c
         maxdia = a(k,k)
         kp1 = k + 1
         maxl = k
c
c        determine the pivot element.
c
         if (k .lt. pl .or. k .ge. pu) go to 190
            do 180 l = kp1, pu
               if (a(l,l) .le. maxdia) go to 170
                  maxdia = a(l,l)
                  maxl = l
  170          continue
  180       continue
  190    continue
c
c        quit if the pivot element is not positive.
c
         if (maxdia .gt. 0.0d0) go to 200
            info = k - 1
c     ......exit
            go to 280
  200    continue
         if (k .eq. maxl) go to 210
c
c           start the pivoting and update piv.
c
            km1 = k - 1
            call dswap(km1,a(1,k),1,a(1,maxl),1)
            a(maxl,maxl) = a(k,k)
            a(k,k) = maxdia
            jp = piv(maxl)
            piv(maxl) = piv(k)
            piv(k) = jp
  210    continue
c
c        reduction step. pivoting is contained across the rows.
c
         work(k) = dsqrt(a(k,k))
         a(k,k) = work(k)
         if (p .lt. kp1) go to 260
         do 250 j = kp1, p
            if (k .eq. maxl) go to 240
               if (j .ge. maxl) go to 220
                  temp = a(k,j)
                  a(k,j) = a(j,maxl)
                  a(j,maxl) = temp
               go to 230
  220          continue
               if (j .eq. maxl) go to 230
                  temp = a(k,j)
                  a(k,j) = a(maxl,j)
                  a(maxl,j) = temp
  230          continue
  240       continue
            a(k,j) = a(k,j)/work(k)
            work(j) = a(k,j)
            temp = -a(k,j)
            call daxpy(j-k,temp,work(kp1),1,a(kp1,j),1)
  250    continue
  260    continue
  270 continue
  280 continue

C Zero the lower triangle
      do i=2,info
        do j=1,i-1
          a(i,j) = 0.0D0
        enddo
      enddo

C decrement the pivots.
      do i=1,p
        piv(i) = piv(i)-1
      enddo
      return
      end

      SUBROUTINE dpotrs_wrap(chol_fac, b, info, n, m, uplo)

cf2py double precision dimension(n,n), intent(in) :: chol_fac
cf2py integer depend(chol_fac), intent(hide)::n=shape(chol_fac,0)
cf2py integer depend(b), intent(hide)::m=shape(b,1)
cf2py optional character intent(in):: uplo='U'
cf2py integer intent(out)::info
cf2py double precision intent(inplace), dimension(n,m)::b
cf2py threadsafe

      DOUBLE PRECISION chol_fac(n,n), b(n,m)
      INTEGER n, info,m

      CHARACTER uplo

      EXTERNAL DPOTRS
! DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO ) Solves triangular system

      call DPOTRS(uplo,n,m,chol_fac,n,b,n,info)

      return
      END

max_stars_repo_name: matthew-brett/pymc

max_issues_repo_name: matthew-brett/pymc

max_forks_repo_name: matthew-brett/pymc


text:
program main
  use iso_c_binding
  implicit none
  integer :: i

  type, bind(c) :: B
    integer(4) :: mb
  end type

  type, bind(c)    :: A
    type(B)        :: tb
    type(B)        :: tb_static_array(5)
    ! generate per array
    type(c_ptr)    :: tb_dyn_array = c_null_ptr
    integer(c_int) :: tb_dyn_array_n1
    integer(c_int) :: tb_dyn_array_lb1
  end type

  interface
     subroutine read_nested_struct(ta) bind(c,name="read_nested_struct")
       import A
       type(A) :: ta
     end subroutine
  end interface

  type(A),target  :: ta
  !
  type(c_ptr)     :: cptr = c_null_ptr
  type(A),pointer :: fptr => null()
  !
  type(B),allocatable,dimension(:) :: tb_dyn_array

  ta%tb%mb = 251

  ! call C routine
  call read_nested_struct(ta)

  ! cast c_ptr to Fortran derived type ptr
  cptr = c_loc(ta)
  call c_f_pointer(cptr,fptr)
  print *, fptr%tb%mb

  ! cast
  allocate(tb_dyn_array(-5:5))

  do i=-5,5
    tb_dyn_array(i)%mb = i
  end do
  !print *, lbound(tb_dyn_array)

  call init_a_tb_dyn_array(ta,tb_dyn_array,lbound(tb_dyn_array,1))
  ! bounds must be passed too. They might get "lost" in subroutine call and reset to 1 (strange)

  call read_nested_struct(ta)
contains

  ! generate
  subroutine init_a_tb_dyn_array(ta,initial,initial_lb1)
     type(A)                     :: ta
     type(B),target,dimension(:) :: initial
     integer                     :: initial_lb1
     !
     !ta%tb_dyn_array_lb1 = lbound(initial,1)
     !print *, lbound(initial,1)
     ta%tb_dyn_array_n1  = size(initial,1)
     ta%tb_dyn_array_lb1 = initial_lb1
     ta%tb_dyn_array     = c_loc(initial)
     !ta%tb_dyn_array     = c_loc(initial(ta%tb_dyn_array_lb1))
  end subroutine

end program

max_stars_repo_name: ROCmSoftwarePlatform/gpufort

max_issues_repo_name: mjklemm/gpufort

max_forks_repo_name: ROCmSoftwarePlatform/gpufort


text:
!! This module contains optimization subroutine and parametric model
module mod_minimize
  !! This module contains optimization subroutine and parametric model
  use mod_constants
  use mod_array
  use mod_optimize
  use mod_optimize_lym

  implicit none

  private

  public :: minimize_spec, minimize

contains

  subroutine minimize_spec(n, m, x, lb, ub, line, dim_v, n_gauss, maxiter, iprint)
    !! Minimize algorithn for a specturm
    implicit none

    integer, intent(in) :: n
    integer, intent(in) :: m
    integer, intent(in) :: dim_v
    integer, intent(in) :: n_gauss, maxiter
    integer, intent(in) :: iprint

    real(xp), intent(in), dimension(:), allocatable :: lb, ub
    real(xp), intent(in), dimension(dim_v) :: line

    real(xp), intent(in), dimension(:), allocatable :: x

    real(xp), parameter    :: factr  = 1.0d+7, pgtol  = 1.0d-5

    character(len=60)      :: task, csave
    logical                :: lsave(4)
    integer                :: isave(44)
    real(xp)               :: f
    real(xp)               :: dsave(29)
    integer,  dimension(:), allocatable  :: nbd, iwa
    real(xp), dimension(:), allocatable  :: g, wa

    real(xp), dimension(:), allocatable  :: residual

    !     Allocate dynamic arrays
    allocate(nbd(n), g(n))
    allocate(iwa(3*n))
    allocate(wa(2*m*n + 5*n + 11*m*m + 8*m))

    allocate(residual(dim_v))
    residual = 0._xp

    ! Init nbd
    nbd = 2

    !     We now define the starting point.
    !     We start the iteration by initializing task.
    task = 'START'

    !     The beginning of the loop
    do while(task(1:2).eq.'FG'.or. task.eq.'NEW_X' .or. task.eq.'START')

       !     This is the call to the L-BFGS-B code.
       call setulb (n, m, x, lb, ub, nbd, f, g, factr, pgtol, wa, iwa, task, iprint, csave, lsave, isave, dsave)

       if (task(1:2) .eq. 'FG') then
          !     Compute function f and gradient g for the sample problem.
          call myresidual(x, line, residual, n_gauss, dim_v)
          f = myfunc_spec(residual)
          call mygrad_spec(n_gauss, g, residual, x, dim_v)

       elseif (task(1:5) .eq. 'NEW_X') then
          !        1) Terminate if the total number of f and g evaluations
          !             exceeds maxiter.
          if (isave(34) .ge. maxiter) &
               task='STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT'

          !        2) Terminate if  |proj g|/(1+|f|) < 1.0d-10.
          if (dsave(13) .le. 1.d-10*(1.0d0 + abs(f))) &
               task='STOP: THE PROJECTED GRADIENT IS SUFFICIENTLY SMALL'
       endif
    !     end of loop do while
    end do

    deallocate(nbd, g)
    deallocate(iwa)
    deallocate(wa)
    deallocate(residual)
  end subroutine minimize_spec

  ! Minimize algorithn for a cube with regularization
  subroutine minimize(n, m, x, lb, ub, cube, n_gauss, dim_v, dim_y, dim_x, lambda_amp, lambda_mu, lambda_sig, &
       lambda_var_amp, lambda_var_mu, lambda_var_sig, lambda_lym_sig, maxiter, kernel, iprint, std_map, lym, c_lym)
    implicit none

    integer, intent(in) :: n
    integer, intent(in) :: m
    integer, intent(in) :: dim_v, dim_y, dim_x
    integer, intent(in) :: n_gauss, maxiter
    integer, intent(in) :: iprint

    real(xp), intent(in) :: c_lym
    real(xp), intent(in) :: lambda_amp, lambda_mu, lambda_sig
    real(xp), intent(in) :: lambda_var_amp, lambda_var_mu, lambda_var_sig
    real(xp), intent(in) :: lambda_lym_sig
    real(xp), intent(in), dimension(:), allocatable :: lb, ub
    real(xp), intent(in), dimension(:,:,:), allocatable :: cube
    real(xp), intent(in), dimension(:,:), allocatable :: kernel
    real(xp), intent(in), dimension(:,:), allocatable :: std_map

    real(xp), intent(in), dimension(:), allocatable :: x

    logical, intent(in) :: lym !! if true --> activate 2-Gaussian decomposition for Lyman alpha nebula emission

    real(xp), parameter    :: factr  = 1.0d+7, pgtol  = 1.0d-5

    character(len=60)      :: task, csave
    logical                :: lsave(4)
    integer                :: isave(44)
    real(xp)               :: f
    real(xp)               :: dsave(29)
    integer,  dimension(:), allocatable  :: nbd, iwa
    real(xp), dimension(:), allocatable  :: g, wa

    !     Allocate dynamic arrays
    allocate(nbd(n), g(n))
    allocate(iwa(3*n))
    allocate(wa(2*m*n + 5*n + 11*m*m + 8*m))

    f = 0._xp
    g = 0._xp

    ! Init nbd
    nbd = 2

    !     We now define the starting point.
    !     We start the iteration by initializing task.
    task = 'START'

    !     The beginning of the loop
    do while(task(1:2).eq.'FG'.or. task.eq.'NEW_X' .or. task.eq.'START')

       !     This is the call to the L-BFGS-B code.
       call setulb (n, m, x, lb, ub, nbd, f, g, factr, pgtol, wa, iwa, task, iprint, csave, lsave, isave, dsave)

       if (task(1:2) .eq. 'FG') then
          !     Compute function f and gradient g for the sample problem.
          if (lym .eqv. .true.) then
             call f_g_cube_fast_lym(f, g, cube, x, dim_v, dim_y, dim_x, n_gauss, kernel, lambda_amp, lambda_mu, lambda_sig, &
                  lambda_var_amp, lambda_var_mu, lambda_var_sig, lambda_lym_sig, std_map, c_lym)
          else
             call f_g_cube_fast(f, g, cube, x, dim_v, dim_y, dim_x, n_gauss, kernel, lambda_amp, lambda_mu, lambda_sig, &
                  lambda_var_amp, lambda_var_mu, lambda_var_sig, std_map)
          end if

       elseif (task(1:5) .eq. 'NEW_X') then
          !        1) Terminate if the total number of f and g evaluations
          !             exceeds maxiter.
          if (isave(34) .ge. maxiter) &
               task='STOP: TOTAL NO. of f AND g EVALUATIONS EXCEEDS LIMIT'

          !        2) Terminate if  |proj g|/(1+|f|) < 1.0d-10.
          if (dsave(13) .le. 1.d-10*(1.0d0 + abs(f))) &
               task='STOP: THE PROJECTED GRADIENT IS SUFFICIENTLY SMALL'
       endif
    !     end of loop do while
    end do

    deallocate(nbd, g)
    deallocate(iwa)
    deallocate(wa)
  end subroutine minimize

end module mod_minimize

max_stars_repo_name: antoinemarchal/ROHSA

max_issues_repo_name: antoinemarchal/ROHSA

max_forks_repo_name: antoinemarchal/ROHSA


text:
!
! what if module name is wrong.
!
      program module
      print 'module'
      end

ext: f


text:
Address(Erma Lane) is a residential street in Royal Oak that is outside of the city limits.

 Intersecting Streets

  Research Park Drive along with Barthel Lane and Morris Lane
  Gay Way
  Harvey Way

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text: I am a Computer Science and Engineering major. I spend ungodly amounts of time in the basement of Kemper Hall, endearingly nicknamed Computer Science Instructional Facility The Dungeon due to working as tech. support for 20 hours a week and also working on my homework projects.


text:
!*****************************************************************************
!
! Rokko: Integrated Interface for libraries of eigenvalue decomposition
!
! Copyright (C) 2012-2019 by Rokko Developers https://github.com/t-sakashita/rokko
!
! Distributed under the Boost Software License, Version 1.0. (See accompanying
! file LICENSE_1_0.txt or copy at http://www.boost.org/license_1_0.txt)
!
!*****************************************************************************

program generate_eigen_vector
  use rokko
  implicit none
  integer :: dim
  type(rokko_eigen_vector) :: vec

  integer :: i

  dim = 10
  print *,"dimension = ", dim

  call rokko_construct(vec, dim)

  do i = 1, dim
     call rokko_set_elem(vec, i, dble(i))
  enddo
  call rokko_print(vec)

  call rokko_destruct(vec)

end program generate_eigen_vector

max_stars_repo_name: t-sakashita/rokko

max_issues_repo_name: t-sakashita/rokko

max_forks_repo_name: t-sakashita/rokko


!    Copyright (C) 2012 The SPEED FOUNDATION
!    Author: Ilario Mazzieri
!
!    This file is part of SPEED.
!
!    SPEED is free software; you can redistribute it and/or modify it
!    under the terms of the GNU Affero General Public License as
!    published by the Free Software Foundation, either version 3 of the
!    License, or (at your option) any later version.
!
!    SPEED is distributed in the hope that it will be useful, but
!    WITHOUT ANY WARRANTY; without even the implied warranty of
!    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!    Affero General Public License for more details.
!
!    You should have received a copy of the GNU Affero General Public License
!    along with SPEED.  If not, see <http://www.gnu.org/licenses/>.

!> @brief Makes Partitioning and writes files *.mpi
!! @author Ilario Mazzieri
!> @date September, 2013
!> @version 1.0

!> @note A new partition is made if elementdomain.mpi is absent in the
!! workig directory or in the folder given in SPEED.input (MPIFILE).

      subroutine MAKE_PARTITION_AND_MPI_FILES()

      use max_var
      use speed_par

      implicit none

      include 'SPEED.MPI'

!***************************************************************************************************************
!     Partitioning

      if (mpi_id.eq.0)  write(*,'(A)') '---------------------Partitioning----------------------'

      if (mpi_id.eq.0) then
              if(len_trim(mpi_file) .ne. 70) then
                 file_part = mpi_file(1:len_trim(mpi_file)) // '/elemdomain.mpi'
              else
                 file_part = 'elemdomain.mpi'
              endif

              inquire(file=file_part,exist=filefound)

!              write(*,*) filefound, file_part, mpi_id
!        read(*,*)

              if(mpi_np .gt. 1) then

           if(filefound .eqv. .TRUE.) then

              if (mpi_id.eq.0)  write(*,'(A)') 'Reading existing partitioning...'
                       unit_part = 400
                 open(unit_part,file=file_part)
                 read(unit_part,*) trash
                 do i = 1, nelem
                       read(unit_part,*) trash, elem_domain(i)
                 enddo
                 close(unit_part)
                 if (mpi_id.eq.0)  write(*,'(A)') 'Read.'

              else

                 write(*,'(A)') 'Making new partitioning...'
                 call MESH_PARTITIONING(mpi_file,nelem,nnod_macro,mpi_np,con,0)
!                call MPI_BARRIER(mpi_comm, mpi_ierr)
                 !sustitute this part adding output elem_domain in mesh partitioning
                 unit_part = 400
                 open(unit_part, file = file_part)
                 read(unit_part,*) trash
                 do i = 1, nelem
                    read(unit_part,*) trash, elem_domain(i)
                 enddo
                 close(unit_part)
                 write(*,'(A)') 'Made.'

              endif
         endif

      endif!if (mpi_id .eq. 0)

      call MPI_BARRIER(mpi_comm, mpi_ierr)
      call MPI_BCAST(elem_domain,nelem,SPEED_INTEGER,0,MPI_COMM_WORLD,mpi_ierr)

!**************************************************************************************************************

      allocate(elem_index(nelem))

      nelem_dom = 0
      do ie = 1,nelem
         elem_index(ie) = 0
         if (elem_domain(ie).eq.mpi_id) then
            nelem_dom = nelem_dom +1
            elem_index(ie) = nelem_dom
         endif
      enddo

!***************************************************************************************************************
!     SPECTRAL CONNECTIVITY --- ONLY FOR MASTER PROCESS ---

      if (mpi_id.eq.0) then
         write(*,'(A)')
         write(*,'(A)')'--------------Making Spectral connectivity-------------'
      endif

      call MPI_BARRIER(mpi_comm, mpi_ierr)

      if(mpi_id .eq. 0) then

        file_mpi = 'nodedomain.mpi';

        if(len_trim(mpi_file) .ne. 70) then
           file_mpi_new = mpi_file(1:len_trim(mpi_file)) // '/' // file_mpi
        else
           file_mpi_new = file_mpi
        endif

        inquire(file=file_mpi_new, exist=filefound)

        if(filefound .eqv. .FALSE.) then

                allocate(node_weight(nnod_macro))

                 call MAKE_WGT_GRID_NODES(nnod_macro,nelem,con,node_weight,nnz_node_weight)

                allocate(node_pointer(0:nnz_node_weight))

                call MAKE_GRID_NODES(nnod_macro,nelem,con,node_weight,nnz_node_weight,node_pointer)

                deallocate(node_weight)

                con_nnz = nelem +1
                do ie = 1,nelem
                     do j = 1,nmat
                        if (tag_mat(j).eq.con(ie,1)) nn = sdeg_mat(j) +1
                     enddo
                     con_nnz = con_nnz + nn*nn*nn +1
                enddo

                allocate(con_spx(0:con_nnz))

                call MAKE_SPX_CON(nelem,con,nmat,tag_mat,sdeg_mat,&
                                      nnz_node_weight,node_pointer,con_nnz,con_spx,nnod)

                deallocate(node_pointer);  allocate(node_weight(nnod))

                call MAKE_WGT_SPX_NODES(nnod,con_nnz,con_spx,node_weight,nnz_node_weight)

                allocate(node_pointer(0:nnz_node_weight))

                call MAKE_SPX_NODES(nnod,con_nnz,con_spx,node_weight,nnz_node_weight,node_pointer)

                deallocate(node_weight)

          write(*,'(A)') 'Made.'

!     SPECTRAL CONNECTIVITY END --- ONLY FOR MASTER PROCESS ---
!***************************************************************************************************************
!***************************************************************************************************************
!     NODE DOMAIN FOR MPI

          allocate(node_domain(nnod))

          do ie = 1, nelem
             do i = con_spx(ie -1) +1, con_spx(ie) -1
                node_domain(con_spx(i)) = elem_domain(ie)
             enddo
          enddo

          file_mpi = 'nodedomain.mpi'; unit_mpi = 40

          if(len_trim(mpi_file) .ne. 70) then
             file_mpi_new = mpi_file(1:len_trim(mpi_file)) // '/' // file_mpi
          else
             file_mpi_new = file_mpi
          endif

          open(unit_mpi,file=file_mpi_new)
          write(unit_mpi,*) nnod
          do i = 1, nnod
             write(unit_mpi,*) node_domain(i)
          enddo
          close(unit_mpi)

          deallocate(node_domain)

!     NODE DOMAIN FOR MPI
!***************************************************************************************************************
!***************************************************************************************************************
!     SPECTRAL LOCAL CONNECTIVITY

          write(*,'(A)')
          write(*,'(A)') '---------------Making local connectivity---------------'

          do ip = 0, mpi_np-1

              nelem_loc = 0
              do in = 1,nelem
                 if(elem_domain(in) .eq. ip) nelem_loc = nelem_loc +1
                    enddo

              con_nnz_loc = nelem_loc +1
              do ie = 1, nelem
                 do j = 1, nmat
                    if ((tag_mat(j).eq. con_spx(con_spx(ie-1))) .and. (elem_domain(ie) .eq. ip)) then
                       nn = sdeg_mat(j) +1
                       con_nnz_loc = con_nnz_loc + nn*nn*nn +1
                    endif
                 enddo
              enddo

              allocate(con_spx_loc(0:con_nnz_loc))

              call MAKE_SPX_CON_LOC(nelem_loc, nelem, elem_domain, &
                        con_nnz,  con_spx, con_nnz_loc, con_spx_loc, &
                        nmat, tag_mat, sdeg_mat, ip)

              file_mpi = 'cons000000.mpi'
              unit_mpi = 40
              if (ip.lt. 10) then
                  write(file_mpi(10:10),'(i1)') ip
              else if (ip .lt. 100) then
                  write(file_mpi(9:10),'(i2)') ip
              else if (ip .lt. 1000) then
                  write(file_mpi(8:10),'(i3)') ip
              else if (ip .lt. 10000) then
                  write(file_mpi(7:10),'(i4)') ip
              else if (ip .lt. 100000) then
                  write(file_mpi(6:10),'(i5)') ip
              else if (ip .lt. 1000000) then
                  write(file_mpi(5:10),'(i6)') ip
              endif

              if(len_trim(mpi_file) .ne. 70) then
                 file_mpi_new = mpi_file(1:len_trim(mpi_file)) // '/' // file_mpi
              else
                 file_mpi_new = file_mpi
              endif

                open(unit_mpi,file=file_mpi_new)
              write(unit_mpi,*) con_nnz_loc
              do i = 0, con_nnz_loc
                       write(unit_mpi,*) con_spx_loc(i)
              enddo
              close(unit_mpi)

              ic = 0
              do in = 1,nelem
                 if (elem_domain(in) .eq. ip) ic = ic +1
              enddo
              nelem_loc = ic

              allocate(local_el_num(nelem_loc))

              ic = 0
              do in = 1,nelem
                 if (elem_domain(in) .eq. ip) then
                    ic = ic +1
                    local_el_num(ic) = in
                 endif
              enddo

                   con_nnz_bc_loc = 0
                   nface_loc = 0

              ic = 0
              if (nface.gt.0) then

                do i = 1,nface

                  call GET_ELEM_FROM_FACE(nnz_node_weight,node_pointer,con_bc(i,2),con_bc(i,3),con_bc(i,5),ie)
                  call GET_INDLOC_FROM_INDGLO(local_el_num, nelem_loc, ie, ic)

                  if (ic .ne. 0) then !I found global index in local_el_num

                    do j = 1,nmat

                      if (tag_mat(j) .eq. con_spx_loc(con_spx_loc(ic-1))) then
                          nn = sdeg_mat(j) +1
                          con_nnz_bc_loc = con_nnz_bc_loc +nn*nn +1
                          nface_loc = nface_loc + 1
                      endif

                    enddo
                  endif

                enddo

                con_nnz_bc_loc = con_nnz_bc_loc + nface_loc + 1

                allocate(con_spx_bc_loc(0:con_nnz_bc_loc))
                con_spx_bc_loc(0) = nface_loc + 1

                call MAKE_SPX_CON_LOC_BOUND(con_nnz_loc,con_spx_loc, nface,con_bc, &
                                         nmat,tag_mat,sdeg_mat,nnz_node_weight,node_pointer, &
                                         con_nnz_bc_loc,con_spx_bc_loc, local_el_num, nelem_loc, &
                                         nface_loc)

              endif

              file_mpi = 'conb000000.mpi'
              unit_mpi = 40
              if (ip.lt. 10) then
                  write(file_mpi(10:10),'(i1)') ip
              else if (ip .lt. 100) then
                  write(file_mpi(9:10),'(i2)') ip
              else if (ip .lt. 1000) then
                  write(file_mpi(8:10),'(i3)') ip
              else if (ip .lt. 10000) then
                  write(file_mpi(7:10),'(i4)') ip
              else if (ip .lt. 100000) then
                  write(file_mpi(6:10),'(i5)') ip
              else if (ip .lt. 1000000) then
                  write(file_mpi(5:10),'(i6)') ip
              endif

              if(len_trim(mpi_file) .ne. 70) then
                 file_mpi_new = mpi_file(1:len_trim(mpi_file)) // '/' // file_mpi
              else
                 file_mpi_new = file_mpi
              endif

              open(unit_mpi,file=file_mpi_new)
              write(unit_mpi,*) con_nnz_bc_loc
              do i = 0, con_nnz_bc_loc
                       write(unit_mpi,*) con_spx_bc_loc(i)
              enddo
              close(unit_mpi)

           deallocate(con_spx_loc, con_spx_bc_loc, local_el_num)

          enddo

          deallocate(con_spx,node_pointer)

        endif

      endif

      if(mpi_id.eq.0) write(*,'(A)') 'Made.'
      call MPI_BARRIER(mpi_comm, mpi_ierr)

!***************************************************************************************************************
!     LOCAL NUMERATION FOR ELEMENTS

      ic = 0
      do in = 1,nelem

         if (elem_domain(in) .eq. mpi_id) then
            ic = ic +1
         endif
      enddo
      nelem_loc = ic

      allocate(local_el_num(nelem_loc))

      ic = 0
      do in = 1,nelem

         if (elem_domain(in) .eq. mpi_id) then
            ic = ic +1
            local_el_num(ic) = in
         endif
      enddo

      call MPI_BARRIER(mpi_comm, mpi_ierr)

!     LOCAL NUMERATION FOR ELEMENTS  END
!***************************************************************************************************************

       file_mpi = 'cons000000.mpi';   unit_mpi = 40

       if (mpi_id .lt. 10) then
          write(file_mpi(10:10),'(i1)') mpi_id
       else if (mpi_id .lt. 100) then
          write(file_mpi(9:10),'(i2)') mpi_id
       else if (mpi_id .lt. 1000) then
          write(file_mpi(8:10),'(i3)') mpi_id
       else if (mpi_id .lt. 10000) then
          write(file_mpi(7:10),'(i4)') mpi_id
       else if (mpi_id .lt. 100000) then
          write(file_mpi(6:10),'(i5)') mpi_id
       else if (mpi_id .lt. 1000000) then
          write(file_mpi(5:10),'(i6)') mpi_id
       endif

       if(len_trim(mpi_file) .ne. 70) then
          file_mpi_new = mpi_file(1:len_trim(mpi_file)) // '/' // file_mpi
       else
          file_mpi_new = file_mpi
       endif

       open(unit_mpi,file=file_mpi_new)
       read(unit_mpi,*) con_nnz_loc
       allocate(con_spx_loc(0:con_nnz_loc))

       do i = 0, con_nnz_loc
            read(unit_mpi,*) con_spx_loc(i)
       enddo
       close(unit_mpi)

       file_mpi = 'conb000000.mpi';  unit_mpi = 40

       if (mpi_id .lt. 10) then
          write(file_mpi(10:10),'(i1)') mpi_id
       else if (mpi_id .lt. 100) then
          write(file_mpi(9:10),'(i2)') mpi_id
       else if (mpi_id .lt. 1000) then
          write(file_mpi(8:10),'(i3)') mpi_id
       else if (mpi_id .lt. 10000) then
          write(file_mpi(7:10),'(i4)') mpi_id
       else if (mpi_id .lt. 100000) then
          write(file_mpi(6:10),'(i5)') mpi_id
       else if (mpi_id .lt. 1000000) then
          write(file_mpi(5:10),'(i6)') mpi_id
       endif

       if(len_trim(mpi_file) .ne. 70) then
          file_mpi_new = mpi_file(1:len_trim(mpi_file)) // '/' // file_mpi
       else
          file_mpi_new = file_mpi
       endif

       open(unit_mpi,file=file_mpi_new)
       read(unit_mpi,*) con_nnz_bc_loc
       allocate(con_spx_bc_loc(0:con_nnz_bc_loc))

      do i = 0, con_nnz_bc_loc
            read(unit_mpi,*) con_spx_bc_loc(i)
      enddo
      close(unit_mpi)

!     SPECTRAL LOCAL CONNECTIVITY
!***************************************************************************************************************
!***************************************************************************************************************
!     COUNTING LOCAL NODES & LOCAL NODE NUMBERING
      if(mpi_id .eq.0) write(*,'(A)')
      if(mpi_id .eq.0) write(*,'(A)') '------------------Counting local nodes-----------------'

      call GET_LOC_NODE_NUM(con_nnz_loc, con_spx_loc, nnod_loc, mpi_id)

      if(mpi_id .eq.0) write(*,'(A)') 'Made.'

      allocate(local_node_num(nnod_loc))
      local_node_num = 0

      if(mpi_id .eq.0) write(*,'(A)')
      if(mpi_id .eq.0) write(*,'(A)') '--------------Make local nodes numbering---------------'

      call MAKE_REN_LOC_NODE(con_nnz_loc, con_spx_loc, nnod_loc, mpi_id, &
                             local_node_num)

!      ic = 0
!!      do im = 1,nmat
!!         nn = sdeg_mat(im) +1
!         do ie = 1,nelem_loc
!            im = con_spx_loc(con_spx_loc(ie -1) +0); nn = sdeg_mat(im) +1
!
!!            if (con_spx_loc(con_spx_loc(ie -1) +0) .eq. tag_mat(im)) then
!                do i = 1, nn*nn*nn
!                   if (find(nnod_loc, local_node_num, con_spx_loc(con_spx_loc(ie -1) + i)) .ne. 1) then
!                      ic = ic + 1
!                      local_node_num(ic) = con_spx_loc(con_spx_loc(ie -1) + i)
!                   endif
!                 enddo
!!             endif
!          enddo
!!       enddo

     if(mpi_id .eq.0) write(*,'(A)') 'Made.'
     if(mpi_id .eq.0) write(*,'(A)')

     file_mpi = 'nodedomain.mpi'; unit_mpi = 40
     if(len_trim(mpi_file) .ne. 70) then
        file_mpi_new = mpi_file(1:len_trim(mpi_file)) // '/' // file_mpi
     else
        file_mpi_new = file_mpi
     endif

     if(mpi_id .eq.0) then
        !start1=MPI_WTIME()
       open(unit_mpi,file=file_mpi_new)
       read(unit_mpi,*) nnod
     endif
     call MPI_BCAST(nnod,1,SPEED_INTEGER,0,MPI_COMM_WORLD,mpi_ierr)

     allocate(node_domain(nnod))
     if(mpi_id .eq.0) then
       do i = 1, nnod
            read(unit_mpi,*) node_domain(i)
       enddo
       close(unit_mpi)
       !write (*,*) "File nodedomain.mpi read, broadcast to slaves"
     endif
     call MPI_BCAST(node_domain,nnod,SPEED_INTEGER,0,MPI_COMM_WORLD,mpi_ierr)

       !write (*,*) "Vector received by process",mpi_id

      ic = 0
      do in = 1, nnod
         if (node_domain(in).eq. mpi_id) ic = ic +1
      enddo
      nnode_dom = ic

  !    write(*,'(A,I3,I8)')'Nodes active on proc : ', mpi_id, ic
      write(*,'(A,I6,I15)') 'Total nodes on Proc. :  ', mpi_id, nnode_dom

      call MPI_BARRIER(mpi_comm, mpi_ierr)

      !if(mpi_id .eq.0) then
         !start2=MPI_WTIME()
         !write (*,*) "Time to send nodedomain.mpi to all &
         !& proc",start2-start1,"sec"
      !endif

!     COUNTING LOCAL NODES & LOCAL NODE NUMBERING
!***********************************************************************************

      end subroutine MAKE_PARTITION_AND_MPI_FILES


text:
      subroutine bigout(bigcrp,iiyear,nowcrp)
c
      include 'pmxcrp.inc'
      include 'pmxelm.inc'
      include 'pmxhil.inc'
      include 'pmxnsl.inc'
      include 'pmxpln.inc'
      include 'pmxpnd.inc'
      include 'pmxres.inc'
      include 'pmxsrg.inc'
      include 'pmxtil.inc'
      include 'pmxtls.inc'
      include 'pntype.inc'
      include 'pxstep.inc'
c
      include 'cavepar.inc'
      include 'cavloss.inc'
      include 'cclim.inc'
      include 'ccons.inc'
      include 'ccover.inc'
      include 'ccrpout.inc'
      include 'ccrpprm.inc'
      include 'ccrpvr1.inc'
      include 'ccrpvr2.inc'
      include 'ccdrain.inc'
c
      include 'cefflen.inc'
c     read: efflen(mxplan)
c
      include 'cends4.inc'
      include 'cenrpa2.inc'
      include 'cffact.inc'

      include 'cflags.inc'
c     read: bigflg
c
      include 'chydrol.inc'
      include 'cirriga.inc'
      include 'cirfurr.inc'
      include 'cparame.inc'
      include 'cparva2.inc'
      include 'cseddet.inc'
      include 'cslpopt.inc'
      include 'csolva2.inc'
      include 'cstruc.inc'
      include 'cupdate.inc'
      include 'cwater.inc'
      include 'cwint.inc'
      include 'ccliyr.inc'
      include 'ccntour.inc'
c
c
c     parameters
c
c
c     local values
c
c
c     static local (hence the `save')
c
      integer tint(8), bigcrp, iiyear, intmin(8), intmax(8),i,kk,nowcrp
      real treal(96), ralmax(96), ralmin(96), watcon, watconf
      save intmin, intmax, ralmin, ralmax
c
      data ralmin /96 * 1e6/
      data ralmax /96 * -1e6/
      data intmin /8 * 1e4/
      data intmax /8 * -1e4/
c
c
c     Daily: UPDATE CURRENT MINIMUMS AND MAXIMUMS AND WRITE
c
      if (bigflg.eq.0) then
c
c       very first time in
c
        if (sdate.eq.1.and.iiyear.eq.1) then
          intmin(1) = 1
          intmax(1) = 1
c
c       increment the date.
c
        else if (iplane.eq.1) then
          intmax(1) = intmax(1) + 1
        end if
c
c
c       determine total soil water content for all the OFEs this date
c
        watcon = 0.0
        watconf = 0.0
        do 10 i = 1, nsl(iplane)
          watcon = watcon + soilw(i,iplane)
          watconf = watconf + soilf(i,iplane)
   10   continue
c
c
c       determine the current minimums and maximums.
c       values which are constant for the hillslope come first.
c
        tint(1) = intmax(1)
c       tint(2) = itype(nowcrp,iplane)
        tint(2) = bigcrp
        do 20 kk = 1, 3
          tint(2+kk) = iresd(kk,iplane)
          tint(5+kk) = iroot(kk,iplane)
   20   continue
c
        treal(1) = prcp * 1000
        treal(2) = avedet
        treal(3) = maxdet
        treal(4) = ptdet
        treal(5) = avedep * (-1.0)
        treal(6) = maxdep * (-1.0)
        treal(7) = ptdep
        treal(8) = avsole
        treal(9) = tmnavg
        treal(10) = tmxavg
        treal(11) = tmin
        treal(12) = tmax
c
        if (irabrv.ne.0) then
          treal(13) = (irdept(iplane) + iraplo(iplane)) * 1000.0
          treal(14) = irapld(iplane) * 1000.0
        else
          treal(13) = 0.0
          treal(14) = 0.0
        end if
c
c  jrf - if contouring is in effect then don't scale runoff because it will not be valid.
c         11/18/2009
c
        if (contrs(nowcrp,iplane).ne.0) then
           treal(15) = runoff(iplane)*1000.
        else
           treal(15) = (runoff(iplane)*efflen(iplane)/totlen(iplane)) *
     1      1000.0
        endif
        treal(16) = irdgdx(iplane)
        treal(17) = canhgt(iplane)
        treal(18) = cancov(iplane)
        treal(19) = lai(iplane)
        treal(20) = inrcov(iplane)
        treal(21) = rilcov(iplane)
        treal(22) = vdmt(iplane)
        treal(23) = rtmass(iplane)
        treal(24) = rtm15(iplane)
        treal(25) = rtm30(iplane)
        treal(26) = rtm60(iplane)
        treal(27) = rtd(iplane)
        treal(28) = rmagt(iplane)
c
        do 30 kk = 1, 3
          treal(28+kk) = rmogt(kk,iplane)
          treal(31+kk) = smrm(kk,iplane)
          treal(34+kk) = rtm(kk,iplane)
   30   continue
c
        treal(38) = avpor(iplane) * 100.0
        treal(39) = avbd(iplane) / 1000.0
        treal(40) = ks(iplane) * 3.6e6
        treal(41) = sm(iplane) * 1000.0
        treal(42) = (es(iplane) + ep(iplane) + eres(iplane)) * 1000.0
c       treal(39) = thetfc(1,iplane)
c       treal(40) = thetdr(1,iplane)
        treal(43) = drainq(iplane)
        treal(44) = solthk(nsl(iplane),iplane) - satdep(iplane)
        treal(45) = effint(iplane) * 3.6e6
        treal(46) = peakro(iplane) * 3.6e6
        treal(47) = effdrn(iplane) / 3600.0
        treal(48) = enrato(iplane)
        treal(49) = (ki(iplane)*kiadjf(iplane)) / 1000000.0
        treal(50) = (kr(iplane)*kradjf(iplane)) * 1000.0
        treal(51) = shcrit(iplane) * tcadjf(iplane)
        treal(52) = width(iplane)
        treal(53) = ep(iplane) * 1000.0
        treal(54) = es(iplane) * 1000.0
        treal(55) = sep(iplane) * 1000.0
        treal(56) = watstr(iplane)
        treal(57) = temstr(iplane)
        treal(58) = watcon * 1000.0
c
        do 40 kk = 1, mxnsl
          treal(58+kk) = soilw(kk,iplane) * 1000.0
   40   continue
c
        treal(69) = rrc(iplane) * 1000.0
        treal(70) = rh(iplane) * 1000.0
        treal(71) = frdp(iplane) * 1000.0
        treal(72) = thdp(iplane) * 1000.0
        treal(73) = snodpy(iplane) * 1000.0
        treal(74) = wmelt(iplane) * 1000.0
        treal(75) = densg(iplane)
        treal(76) = frccov(iplane)
        treal(77) = frlive(iplane)
        treal(78) = frctrl(iplane)
        treal(79) = frcteq(iplane)
        treal(80) = frrres(iplane)
        treal(81) = fribas(iplane)
        treal(82) = frican(iplane)
        treal(83) = daydis(iplane)
        treal(84) = ofelod(iplane)
        treal(85) = eres(iplane)*1000.0

        treal(86) = watconf * 1000.0
        do 45 kk = 1, mxnsl
          treal(86+kk) = soilf(kk,iplane) * 1000.0
   45   continue

        effdrn(iplane) = 0.0
        temstr(iplane) = 1.0
        frccov(iplane) = 0.0
        frlive(iplane) = 0.0
        frctrl(iplane) = 0.0
        frcteq(iplane) = 0.0
        frrres(iplane) = 0.0
        fribas(iplane) = 0.0
        frican(iplane) = 0.0
c       ofelod(iplane) = 0.0
c
c
c       determine the minimums and maximums.
c
        do 50 i = 1, 96
          call mxreal(treal(i),ralmin(i),ralmax(i))
   50   continue
        do 60 i = 1, 8
          call mxint(tint(i),intmin(i),intmax(i))
   60   continue
c
c
c       write the daily information
c
        write (40,1000) tint(1), (treal(i),i=1,83), (tint(i),i = 2,8),
     1         (treal(i),i=84,96)
c
c
c     WEPP has completed - append min/max values.
c
      else
c
        write (40,*) '#'
        write (40,*) '#       Minimum/Maximum values:'
        write (40,*) '#'
c
        write (40,1000)intmin(1),(ralmin(i),i=1,83),(intmin(i),i = 2,8),
     1        (ralmin(i),i=84,96)
        write (40,1000)intmax(1),(ralmax(i),i=1,83),(intmax(i),i = 2,8),
     1        (ralmax(i),i=84,96)
c
      end if
      return
c
 1000 format (i6,83(1x,f10.5),7(1x,i2),13(1x,f10.5))
      end

ext: for

max_stars_repo_name: jarad/dep

max_issues_repo_name: jarad/dep

max_forks_repo_name: akrherz/idep


text:
MODULE common
!==============================================================================
!
! MODULE: letkf_local
!
! PURPOSE: General constants and procedures
!
! USES:
!   This module calls 'SFMT.f90' for randomization functions
!
! DESCRIPTION:
!   This is the main program for the letkf data assimilation.
!
! REVISION HISTORY:
!   07/20/2004 Takemasa Miyoshi  created
!   01/23/2009 Takemasa Miyoshi  modified for SFMT
!   04/26/2011 Steve Penny converted to OCEAN for use with MOM4
!   07/22/2015 Steve Penny removed unused subroutines
!
!-------------------------------------------------------------------------------
! $Authors: Steve Penny, Takemasa Miyoshi $
!===============================================================================

  IMPLICIT NONE
  PUBLIC

!-----------------------------------------------------------------------
! Variable size definitions
!-----------------------------------------------------------------------
  INTEGER,PARAMETER :: r_size=kind(0.0d0)
  INTEGER,PARAMETER :: r_dble=kind(0.0d0)
  INTEGER,PARAMETER :: r_sngl=kind(0.0e0)
  INTEGER,PARAMETER :: slen=512

!-----------------------------------------------------------------------
! Constants
!-----------------------------------------------------------------------
  REAL(r_size),PARAMETER :: pi=3.1415926535d0
  REAL(r_size),PARAMETER :: gg=9.81d0
  REAL(r_size),PARAMETER :: rd=287.0d0
  REAL(r_size),PARAMETER :: cp=1005.7d0
  REAL(r_size),PARAMETER :: re=6371.3d3      !STEVE: i.e. d3 => re is in meters
  REAL(r_size),PARAMETER :: r_omega=7.292d-5
  REAL(r_size),PARAMETER :: t0c=273.15d0
  REAL(r_size),PARAMETER :: undef=-9.99d33

!-----------------------------------------------------------------------
! STEVE:
! LETKF Parameters
!-----------------------------------------------------------------------
  REAL(r_size),PARAMETER :: sigma_b = 0.1d0 !1.0d0 !0.1d0

CONTAINS

SUBROUTINE com_filter_lanczos(ndim,fc,var)
!-----------------------------------------------------------------------
! Lanczos Filter (Low-pass) with cyclic boundary
!-----------------------------------------------------------------------
  IMPLICIT NONE

  INTEGER,INTENT(IN) :: ndim
  REAL(r_size),INTENT(IN) :: fc    ! critical frequency in [0,pi]
  REAL(r_size),INTENT(INOUT) :: var(ndim)

  INTEGER,PARAMETER :: lresol=10

  REAL(r_size) :: weight(-lresol:lresol)
  REAL(r_size) :: varwk(1-lresol:ndim+lresol)
  REAL(r_size) :: rl,rlresol
  INTEGER :: i,l
!
! Weight
!
  rlresol = REAL(lresol,r_size)
  DO l=-lresol,-1
    rl = REAL(l,r_size)
    weight(l) = SIN(fc*rl) * SIN(pi*rl/rlresol) &
      & * rlresol / pi / rl / pi / rl
  END DO
  DO l=1,lresol
    rl = REAL(l,r_size)
    weight(l) = SIN(fc*rl) * SIN(pi*rl/rlresol) &
      & * rlresol / pi / rl / pi / rl
  END DO
  weight(0) = fc / pi
!
! Cyclic boundary
!
  DO i=0,1-lresol,-1
    varwk(i) = var(ndim+i)
  END DO
  DO i=ndim+1,ndim+lresol
    varwk(i) = var(i-ndim)
  END DO
  varwk(1:ndim) = var(1:ndim)
!
! Filter
!
  var = 0.0d0
  DO i=1,ndim
    DO l=-lresol,lresol
      var(i) = var(i) + weight(l) * varwk(i+l)
    END DO
  END DO

END SUBROUTINE com_filter_lanczos

SUBROUTINE com_rand(ndim,var)
!-----------------------------------------------------------------------
! RAND (random number with uniform distribution)
!-----------------------------------------------------------------------
  IMPLICIT NONE

  INTEGER,INTENT(IN) :: ndim
  REAL(r_size),INTENT(OUT) :: var(1:ndim)
  REAL(r_dble) :: genrand_res53
  INTEGER :: idate(8)
  INTEGER :: i,iseed
  LOGICAL,SAVE :: first=.true.

  IF (first) THEN
    CALL DATE_AND_TIME(VALUES=idate)
    iseed = idate(8) + idate(7)*1000
    CALL init_gen_rand(iseed)
    first=.false.
  END IF

  DO i=1,ndim
    var(i) = genrand_res53()
  END DO

END SUBROUTINE com_rand

SUBROUTINE com_randn(ndim,var)
!-----------------------------------------------------------------------
! RANDN (random number with normal distribution)
!-----------------------------------------------------------------------
  IMPLICIT NONE

  INTEGER,INTENT(IN) :: ndim
  REAL(r_size),INTENT(OUT) :: var(1:ndim)
  REAL(r_size) :: rnd(2)
  REAL(r_dble) :: genrand_res53
  INTEGER :: idate(8)
  INTEGER :: i,iseed
  LOGICAL,SAVE :: first=.true.

  IF (first) THEN
    CALL DATE_AND_TIME(VALUES=idate)
    iseed = idate(8) + idate(7)*1000
    CALL init_gen_rand(iseed)
    first=.false.
  END IF

  IF( MOD(ndim,2)==0 ) THEN
    DO i=1,ndim/2
      rnd(1) = genrand_res53()
      rnd(2) = genrand_res53()
      var(i*2-1) = sqrt( -2.0d0 * log( rnd(1) ) ) * sin( 2.0d0*pi*rnd(2) )
      var(i*2) = sqrt( -2.0d0 * log( rnd(1) ) ) * cos( 2.0d0*pi*rnd(2) )
    END DO
  ELSE
    DO i=1,(ndim-1)/2
      rnd(1) = genrand_res53()
      rnd(2) = genrand_res53()
      var(i*2-1) = sqrt( -2.0d0 * log( rnd(1) ) ) * sin( 2.0d0*pi*rnd(2) )
      var(i*2) = sqrt( -2.0d0 * log( rnd(1) ) ) * cos( 2.0d0*pi*rnd(2) )
    END DO
    rnd(1) = genrand_res53()
    rnd(2) = genrand_res53()
    var(ndim) = sqrt( -2.0d0 * log( rnd(1) ) ) * sin( 2.0d0*pi*rnd(2) )
  END IF

END SUBROUTINE com_randn

SUBROUTINE com_distll(ndim,alon,alat,blon,blat,dist)
!-----------------------------------------------------------------------
! DISTANCE BETWEEN TWO POINTS (LONa,LATa)-(LONb,LATb)
!-----------------------------------------------------------------------
  IMPLICIT NONE
  INTEGER,INTENT(IN) :: ndim
  REAL(r_size),INTENT(IN) :: alon(ndim)
  REAL(r_size),INTENT(IN) :: alat(ndim)
  REAL(r_size),INTENT(IN) :: blon(ndim)
  REAL(r_size),INTENT(IN) :: blat(ndim)
  REAL(r_size),INTENT(OUT) :: dist(ndim)
  REAL(r_size),PARAMETER :: r180=1.0d0/180.0d0
  REAL(r_size) :: lon1,lon2,lat1,lat2
  REAL(r_size) :: cosd(ndim)
  INTEGER :: i

  DO i=1,ndim
    lon1 = alon(i) * pi * r180
    lon2 = blon(i) * pi * r180
    lat1 = alat(i) * pi * r180
    lat2 = blat(i) * pi * r180

    cosd(i) = SIN(lat1)*SIN(lat2) + COS(lat1)*COS(lat2)*COS(lon2-lon1)
    cosd(i) = MIN( 1.d0,cosd(i))
    cosd(i) = MAX(-1.d0,cosd(i))

    dist(i) = ACOS( cosd(i) ) * re
  END DO

END SUBROUTINE com_distll

SUBROUTINE com_distll_1(alon,alat,blon,blat,dist)
!-----------------------------------------------------------------------
! DISTANCE BETWEEN TWO POINTS (LONa,LATa)-(LONb,LATb)
!-----------------------------------------------------------------------
  IMPLICIT NONE
  REAL(r_size),INTENT(IN) :: alon
  REAL(r_size),INTENT(IN) :: alat
  REAL(r_size),INTENT(IN) :: blon
  REAL(r_size),INTENT(IN) :: blat
  REAL(r_size),INTENT(OUT) :: dist
  REAL(r_size),PARAMETER :: r180=1.0d0/180.0d0
  REAL(r_size) :: lon1,lon2,lat1,lat2
  REAL(r_size) :: cosd

  lon1 = alon * pi * r180
  lon2 = blon * pi * r180
  lat1 = alat * pi * r180
  lat2 = blat * pi * r180

  cosd = SIN(lat1)*SIN(lat2) + COS(lat1)*COS(lat2)*COS(lon2-lon1)
  cosd = MIN( 1.d0,cosd)
  cosd = MAX(-1.d0,cosd)

  dist = ACOS( cosd ) * re

END SUBROUTINE com_distll_1

SUBROUTINE com_interp_spline(ndim,x,y,n,x5,y5)
!-----------------------------------------------------------------------
! Cubic spline interpolation
!   [Reference:] Akima, H., 1970: J. ACM, 17, 589-602.
!-----------------------------------------------------------------------
  IMPLICIT NONE
  INTEGER,INTENT(IN) :: ndim         ! number of grid points
  REAL(r_size),INTENT(IN) :: x(ndim) ! coordinate
  REAL(r_size),INTENT(IN) :: y(ndim) ! variable
  INTEGER,INTENT(IN) :: n            ! number of targets
  REAL(r_size),INTENT(IN) :: x5(n)   ! target coordinates
  REAL(r_size),INTENT(OUT) :: y5(n)  ! target values
  INTEGER :: i,j,m
  REAL(r_size) :: dydx(5),ddydx(4),t(2),dx21,dx
  REAL(r_size) :: wk

  TGT: DO j=1,n
    DO i=1,ndim
      IF(x5(j) == x(i)) THEN
        y5(j) = y(i)
        CYCLE TGT
      END IF
      IF(x5(j) < x(i)) EXIT
    END DO
!       i-3   i-2   i-1    i    i+1   i+2
!     ---+-----+-----+---*-+-----+-----+---
!dydx       1     2     3     4     5
!ddydx         1     2     3     4
!t                   1     2
    IF(i==2) THEN
      DO m=3,5
        dydx(m) = (y(i-3+m)-y(i-4+m)) / (x(i-3+m)-x(i-4+m))
      END DO
      dydx(2) = 2.0d0*dydx(3) - dydx(4)
      dydx(1) = 2.0d0*dydx(2) - dydx(3)
    ELSE IF(i==3) THEN
      DO m=2,5
        dydx(m) = (y(i-3+m)-y(i-4+m)) / (x(i-3+m)-x(i-4+m))
      END DO
      dydx(1) = 2.0d0*dydx(2) - dydx(3)
    ELSE IF(i==ndim) THEN
      DO m=1,3
        dydx(m) = (y(i-3+m)-y(i-4+m)) / (x(i-3+m)-x(i-4+m))
      END DO
      dydx(4) = 2.0d0*dydx(3) - dydx(2)
      dydx(5) = 2.0d0*dydx(4) - dydx(3)
    ELSE IF(i==ndim-1) THEN
      DO m=1,4
        dydx(m) = (y(i-3+m)-y(i-4+m)) / (x(i-3+m)-x(i-4+m))
      END DO
      dydx(5) = 2.0d0*dydx(4) - dydx(3)
    ELSE
      DO m=1,5
        dydx(m) = (y(i-3+m)-y(i-4+m)) / (x(i-3+m)-x(i-4+m))
      END DO
    END IF
    DO m=1,4
      ddydx(m) = ABS(dydx(m+1) - dydx(m))
    END DO
    DO m=1,2
      wk = ddydx(m+2) + ddydx(m)
      IF(wk == 0) THEN
        t(m) = 0.0d0
      ELSE
        t(m) = (ddydx(m+2)*dydx(m+1)+ddydx(m)*dydx(m+2))/wk
      END IF
    END DO
    dx21 = x(i)-x(i-1)
    dx = x5(j) - x(i-1)
    y5(j) = y(i-1) &
        & + dx*t(1) &
        & + dx*dx*(3.0d0*dydx(3)-2.0d0*t(1)-t(2))/dx21 &
        & + dx*dx*dx*(t(1)+t(2)-2.0d0*dydx(3))/dx21/dx21
  END DO TGT

END SUBROUTINE com_interp_spline

END MODULE common

max_stars_repo_name: GEOS-ESM/Ocean-LETKF

max_issues_repo_name: GEOS-ESM/Ocean-LETKF

max_forks_repo_name: GEOS-ESM/Ocean-LETKF


text:
!-----------------------------------------------------------------------
!     Subroutine:    stellopt_prof_to_vmec(iflag)
!     Authors:       S. Lazerson (lazerson@pppl.gov)
!     Date:          11/20/2012
!     Description:   This subroutine takes various profile arrays and
!                    constructs the am_aux and ac_aux arrays.  It begins
!                    by redefining the s vector over 25 points
!                    equidistant in toroidal flux.  Then it splines
!                    each quantity to the new grid before recomposing
!                    the quantities onto the new s mesh.
!                    It has also been modified to setup the AH and AT
!                    arrays in order to properly define toroidal
!                    rotation.
!-----------------------------------------------------------------------
      SUBROUTINE stellopt_prof_to_vmec(file_str,iflag)
!-----------------------------------------------------------------------
!     Libraries
!-----------------------------------------------------------------------
      USE stel_kinds, ONLY: rprec
      USE stellopt_vars
      USE stellopt_targets
      USE equil_utils
      USE vmec_input, ONLY: am_aux_s, am_aux_f, ac_aux_s, ac_aux_f, &
                            ah_aux_s, ah_aux_f, at_aux_s, at_aux_f, &
                            pmass_type, pcurr_type, ph_type, pt_type
      USE safe_open_mod
      IMPLICIT NONE

!-----------------------------------------------------------------------
!     Input Variables
!        iflag   Error flag
!----------------------------------------------------------------------
!      CHARACTER(256),INTENT(in) :: file_str   !Changed due to bounds check issue in init_stellopt
      CHARACTER(*),INTENT(in) :: file_str
      INTEGER, INTENT(inout)   :: iflag
!-----------------------------------------------------------------------
!     Local Variables
!        ier         Error flag
!        iunit       File unit number
!----------------------------------------------------------------------
      LOGICAL ::  lnew_am, lnew_ac, lnew_at, lnew_ah, lno_file, lnew_diag
      INTEGER ::  dex_ne, dex_te, dex_ti, dex_zeff, &
                  dex_beamj, dex_bootj, dex_am, dex_ac, &
                  dex_ah, dex_at
      INTEGER ::  ier,ik, iunit, dex
      REAL(rprec) :: emis_xics_temp, w3_xics_temp, phi_temp
      REAL(rprec), ALLOCATABLE :: ne_temp(:), zeff_temp(:), te_temp(:), ti_temp(:),&
                                  th_temp(:), beamj_temp(:), bootj_temp(:), dump_temp(:)

      INTEGER, PARAMETER :: new_prof=50

!----------------------------------------------------------------------
!     BEGIN SUBROUTINE
!----------------------------------------------------------------------
      lno_file = .false.
      IF (iflag==-327) lno_file = .true.
      iflag = 0

      ! Setup Splines if necessary (these also free)
      dex = MINLOC(phi_aux_s(2:),DIM=1)
      IF (dex > 2) CALL setup_prof_spline(phi_spl,dex,phi_aux_s(1:dex),phi_aux_f(1:dex),ier)
      dex = MINLOC(ne_aux_s(2:),DIM=1)
      IF (dex > 2) CALL setup_prof_spline(ne_spl,dex,ne_aux_s(1:dex),ne_aux_f(1:dex),ier)
      dex = MINLOC(te_aux_s(2:),DIM=1)
      IF (dex > 2) CALL setup_prof_spline(te_spl,dex,te_aux_s(1:dex),te_aux_f(1:dex),ier)
      dex = MINLOC(ti_aux_s(2:),DIM=1)
      IF (dex > 2) CALL setup_prof_spline(ti_spl,dex,ti_aux_s(1:dex),ti_aux_f(1:dex),ier)
      dex = MINLOC(th_aux_s(2:),DIM=1)
      IF (dex > 2) CALL setup_prof_spline(th_spl,dex,th_aux_s(1:dex),th_aux_f(1:dex),ier)
      dex = MINLOC(nustar_s(2:),DIM=1)
      IF (dex > 2) CALL setup_prof_spline(nustar_spl,dex,nustar_s(1:dex),nustar_f(1:dex),ier)
      dex = MINLOC(zeff_aux_s(2:),DIM=1)
      IF (dex > 2) CALL setup_prof_spline(zeff_spl,dex,zeff_aux_s(1:dex),zeff_aux_f(1:dex),ier)
      dex = MINLOC(ah_aux_s(2:),DIM=1)
      IF (dex > 2) CALL setup_prof_spline(ah_spl,dex,ah_aux_s(1:dex),ah_aux_f(1:dex),ier)
      dex = MINLOC(emis_xics_s(2:),DIM=1)
      IF (dex > 2) CALL setup_prof_spline(emis_xics_spl,dex,emis_xics_s(1:dex),emis_xics_f(1:dex),ier)
!      dex = MINLOC(omega_spl_s(2:),DIM=1) ! this is a vmec varaible
!      IF (dex > 4) CALL setup_prof_spline(omega_spl,dex,rho(1:dex),omega(1:dex),ier)

      ! First get indexes of component arrays
      dex_ne    = MINLOC(ne_aux_s(2:),DIM=1)
      dex_te    = MINLOC(te_aux_s(2:),DIM=1)
      dex_ti    = MINLOC(ti_aux_s(2:),DIM=1)
      dex_zeff  = MINLOC(zeff_aux_s(2:),DIM=1)
      dex_beamj = MINLOC(beamj_aux_s(2:),DIM=1)
      dex_bootj = MINLOC(bootj_aux_s(2:),DIM=1)
      dex_ah    = MINLOC(ah_aux_s(2:),DIM=1)
      dex_at    = MINLOC(at_aux_s(2:),DIM=1)
      ! Now check if we need new am and ac arrays
      lnew_am = .false.; lnew_ac = .false.; lnew_diag = .false.
      lnew_ah = .false.; lnew_at = .false.
      IF (ANY(ne_opt /= 0.0_rprec)) lnew_am = .true.
      IF (ANY(te_opt /= 0.0_rprec)) lnew_am = .true.
      IF (ANY(ti_opt /= 0.0_rprec)) lnew_am = .true.
      IF (ANY(th_opt /= 0.0_rprec)) lnew_am = .true.
      IF (ANY(bootj_aux_f /= 0.0_rprec)) lnew_ac = .true.
      IF (ANY(beamj_aux_f /= 0.0_rprec)) lnew_ac = .true.
      IF (ANY(emis_xics_f /= 0.0_rprec)) lnew_diag = .true.
      IF (ANY(phi_aux_f /= 0.0_rprec))   lnew_diag = .true.
      IF (dex_te > 3) lnew_am = .true.
      IF (dex_ti > 3) lnew_am = .true.
      !IF (dex_beamj > 1) lnew_ac = .true.
      !IF (dex_bootj > 1) lnew_ac = .true.
      IF (TRIM(equil_type)=='flow' .or. TRIM(equil_type)=='satire') THEN
         lnew_ah = .true.
         lnew_at = .true.
      END IF

      ! Now create the new arrays
      IF (lnew_am) THEN
         iunit = 66
         IF (.not. lno_file) CALL safe_open(iunit,iflag,TRIM('tprof.'//TRIM(file_str)),'unknown','formatted')
         !WRITE(iunit,*) am_aux_s
         am_aux_s(:) = -1
         DO ik = 1, new_prof
            am_aux_s(ik) = REAL(ik-1)/REAL(new_prof-1)
         END DO
         am_aux_s(1) = 0.0_rprec
         am_aux_s(new_prof) = 1.0_rprec
         dex_am = MINLOC(am_aux_s(2:),DIM=1)
         !WRITE(iunit,*) am_aux_s
         ! NE
         ALLOCATE(ne_temp(1:dex_am))
         ne_temp = 1.0_rprec
         DO ik = 1, dex_am
            ier = 0
            CALL get_equil_ne(am_aux_s(ik),TRIM(ne_type),ne_temp(ik),ier)
            IF (ier /= 0) ne_temp(ik) = 1.0D18
         END DO
         ier = 0
         ! TE
         ALLOCATE(te_temp(1:dex_am))
         te_temp = 0.0_rprec
         DO ik = 1, dex_am
            ier = 0
            CALL get_equil_te(am_aux_s(ik),TRIM(te_type),te_temp(ik),ier)
         END DO
         ! TI
         ALLOCATE(ti_temp(1:dex_am))
         ti_temp = 0.0_rprec
         DO ik = 1, dex_am
            ier = 0
            CALL get_equil_ti(am_aux_s(ik),TRIM(ti_type),ti_temp(ik),ier)
         END DO
         ! ZEFF
         ALLOCATE(zeff_temp(1:dex_am))
         zeff_temp = 1.0_rprec
         DO ik = 1, dex_am
            ier = 0
            CALL get_equil_zeff(am_aux_s(ik),TRIM(zeff_type),zeff_temp(ik),ier)
            zeff_temp(ik) = MAX(zeff_temp(ik),1.0)
         END DO
         pmass_type = 'akima_spline'
         ! Create p
         DO ik = 1, dex_am
            am_aux_f(ik) = ne_temp(ik) * ec * (te_temp(ik) + ti_temp(ik)/zeff_temp(ik))
            IF (am_aux_f(ik) < 0.0_rprec) am_aux_f(ik) = 0.0_rprec
         END DO
         !WRITE(iunit,*) am_aux_s
         !WRITE(iunit,*) am_aux_f
         ! Now make table
         IF (.not. lno_file) THEN
            WRITE(iunit,*) 's     ne     te     ti     zeff     p'
            DO ik = 1, dex_am
               WRITE(iunit,'(7es12.4)') am_aux_s(ik),ne_temp(ik),te_temp(ik),ti_temp(ik),zeff_temp(ik),am_aux_f(ik)
            END DO
            CLOSE(iunit)
         END IF
      END IF
      !  SINCE ROTATION does not depend on anything I don't think we need to do
      !  this.  -SAL 12/2/13
      !IF (lnew_ah) THEN
      !   iunit = 66
      !   CALL safe_open(iunit,iflag,TRIM('rotprof.'//TRIM(proc_string)),'unknown','formatted')
      !   ah_aux_s(:) = -1
      !   DO ik = 1, new_prof
      !      ah_aux_s(ik) = REAL(ik-1)/REAL(new_prof-1)
      !   END DO
      !   ah_aux_s(1) = 0.0
      !   ah_aux_s(new_prof) = 1.0
      !   dex_ah = MINLOC(ah_aux_s(2:),DIM=1)
      !   ph_type = 'akima_spline'
      !   DO ik = 1, dex_ah
      !      ier = 0
      !      CALL get_equil_omega(ah_aux_s(ik),ah_aux_f(ik),ier)
      !   END DO
      !   ! Now make table
      !   WRITE(iunit,*) 's     rot'
      !   DO ik = 1, dex_am
      !      WRITE(iunit,'(7es12.4)') ah_aux_s(ik),ah_aux_f(ik)
      !   END DO
      !   CLOSE(iunit)
      !END IF
      IF (lnew_at) THEN
         !iunit = 66
         !CALL safe_open(iunit,iflag,TRIM('rotprof.'//TRIM(proc_string)),'unknown','formatted')
         at_aux_s(:) = -1
         DO ik = 1, new_prof
            at_aux_s(ik) = REAL(ik-1)/REAL(new_prof-1)
         END DO
         at_aux_s(1) = 0.0_rprec
         at_aux_s(new_prof) = 1.0_rprec
         dex_at = MINLOC(at_aux_s(2:),DIM=1)
         pt_type = 'akima_spline'
         DO ik = 1, dex_at
            ier = 0
            CALL get_equil_ti(at_aux_s(ik),TRIM(ti_type),at_aux_f(ik),ier)
         END DO
         ! Now make table
         !WRITE(iunit,*) 's     rot'
         !DO ik = 1, dex_am
         !   WRITE(iunit,'(7es12.4)') ah_aux_s(ik),ah_aux_f(ik)
         !END DO
         !CLOSE(iunit)
      END IF

      IF (lnew_ac) THEN
         iunit = 66
         IF (.not. lno_file) CALL safe_open(iunit,iflag,TRIM('jprof.'//TRIM(file_str)),'unknown','formatted')
         !WRITE(iunit,*) ac_aux_s
         ac_aux_s(:) = -1
         DO ik = 1, new_prof
            ac_aux_s(ik) = REAL(ik-1)/REAL(new_prof-1)
         END DO
         ac_aux_s(1) = 0.0_rprec
         ac_aux_s(new_prof) = 1.0_rprec
         dex_ac = MINLOC(ac_aux_s(2:),DIM=1)
         !WRITE(iunit,*) ac_aux_s
         ! BEAMJ
         ALLOCATE(beamj_temp(1:dex_ac))
         beamj_temp = 0.0_rprec
         DO ik = 1, dex_ac
            ier = 0
            CALL get_equil_beamj(ac_aux_s(ik),beamj_temp(ik),ier)
         END DO
         ! BEAMJ
         ALLOCATE(bootj_temp(1:dex_ac))
         bootj_temp = 0.0_rprec
         DO ik = 1, dex_ac
            ier = 0
            CALL get_equil_bootj(ac_aux_s(ik),bootj_temp(ik),ier)
         END DO
         pcurr_type = 'akima_spline_ip'
         ! Create j
         DO ik = 1, dex_ac
            ac_aux_f(ik) = beamj_temp(ik) + bootj_temp(ik)
         END DO
         !WRITE(iunit,*) ac_aux_s
         !WRITE(iunit,*) ac_aux_f
         IF (.not. lno_file) THEN
            WRITE(iunit,*) 's    j_beam     j_bs     j_tot'
            DO ik = 1, dex_ac
               WRITE(iunit,'(4es12.4)') ac_aux_s(ik),beamj_temp(ik),bootj_temp(ik),ac_aux_f(ik)
            END DO
            CLOSE(iunit)
         END IF
      END IF

      ! Free Arrays
      IF (ALLOCATED(ne_temp)) DEALLOCATE(ne_temp)
      IF (ALLOCATED(zeff_temp)) DEALLOCATE(zeff_temp)
      IF (ALLOCATED(te_temp)) DEALLOCATE(te_temp)
      IF (ALLOCATED(ti_temp)) DEALLOCATE(ti_temp)
      IF (ALLOCATED(bootj_temp)) DEALLOCATE(bootj_temp)
      IF (ALLOCATED(beamj_temp)) DEALLOCATE(beamj_temp)

      ! Output a diagnostic table
      ! Now make table
      IF (lnew_diag) THEN
         iunit = 68
         ALLOCATE(ne_temp(new_prof)) ! Use ne_temp for s
         ne_temp = -1
         DO ik = 1, new_prof
            ne_temp(ik) = REAL(ik-1)/REAL(new_prof-1)
         END DO
         ne_temp(1) = 0.0_rprec
         ne_temp(new_prof) = 1.0_rprec
         IF (.not. lno_file) THEN
            CALL safe_open(iunit,iflag,TRIM('dprof.'//TRIM(file_str)),'unknown','formatted')
            WRITE(iunit,*) 's     emis_xics       phi'
            DO ik = 1, new_prof
               CALL get_equil_emis_xics(ne_temp(ik),TRIM(emis_xics_type),emis_xics_temp,ier)
               CALL get_equil_phi(ne_temp(ik),TRIM(phi_type),phi_temp,ier)
               WRITE(iunit,'(3es12.3)') ne_temp(ik),emis_xics_temp,phi_temp
            END DO
            CLOSE(iunit)
         END IF
         DEALLOCATE(ne_temp)
      END IF

      RETURN
!----------------------------------------------------------------------
!     END SUBROUTINE
!----------------------------------------------------------------------
      END SUBROUTINE stellopt_prof_to_vmec

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
The following businesses have high chairs for Parenting parents and caregivers with young children.

  Burgers and Brew
  Cafe Bernardo
  Crepeville
  Dos Coyotes
  El Mariachi
  Hoa Viet (kind of scary)
  Paesanos
  Sophias Thai Kitchen
  Sushi Unlimited
  Taqueria Guadalajara

But these businesses do NOT have any.

  Mishkas Cafe
  Redrum Burger (they do have one booster seat)
  Swirl

Note: Some places have high chairs, but look like theyre going to fall apart, or are missing the middle strap (goes between the legs), or have broken or missing belt buckles.

See also: Changing Tables


text:
*
* METRo : Model of the Environment and Temperature of Roads
* METRo is Free and is proudly provided by the Government of Canada
* Copyright (C) Her Majesty The Queen in Right of Canada, Environment Canada, 2006

*  Questions or bugs report: metro@ec.gc.ca
*  METRo repository: https://gna.org/projects/metro/
*  Documentation: http://documentation.wikia.com/wiki/METRo
*
*
* Code contributed by:
*  Yves Delage - Recherche en prevision numerique (RPN)
*  Miguel Tremblay - Canadian meteorological center
*
*  $LastChangedDate$
*  $LastChangedRevision$
*
************************************************************************
*  This program is free software; you can redistribute it and/or modify
*  it under the terms of the GNU General Public License as published by
*  the Free Software Foundation; either version 2 of the License, or
*  (at your option) any later version.
*
*  This program is distributed in the hope that it will be useful,
*  but WITHOUT ANY WARRANTY; without even the implied warranty of
*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
*  GNU General Public License for more details.
*
*  You should have received a copy of the GNU General Public License
*  along with this program; if not, write to the Free Software
*  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*
*
*
************************************************************************

***S/P  FLXSURFZ
*
      SUBROUTINE FLXSURFZ(CMU, CTU, RIB, FTEMP, FVAP, ILMO,
     X                    UE, FCOR, TA , QA , ZU, ZT, VA,
     Y                    TG , QG , H , Z0 , Z0T,
     %                    LZZ0, LZZ0T, FM, FH,N, ITER )
      IMPLICIT NONE
      INTEGER N,ITER(N)
      REAL CMU(N),CTU(N),RIB(N),ILMO(N)
      REAL FTEMP(N),FVAP(N)
      DOUBLE PRECISION FCOR(N)
      DOUBLE PRECISION TA(N),QA(N),ZU(N),ZT(N),VA(N)
      DOUBLE PRECISION Z0(N), Z0T(N)
      REAL TG(N),QG(N),H(N),UE(N)
      REAL LZZ0(N),LZZ0T(N)
      REAL fm(N),fh(N)
*
*Author
*          Y.Delage (Jul 1990)
*Revision
* 001      G. Pellerin (Jun 94) New function for unstable case
* 002      G. Pellerin (Jui 94) New formulation for stable case
* 003      B. Bilodeau (Nov 95) Replace VK by KARMAN
* 004      M. Desgagne (Dec 95) Add safety code in function ff
*                               and ensures that RIB is non zero
* 005      R. Sarrazin (Jan 96) Correction for H
* 006      C. Girard (Nov 95) - Diffuse T instead of Tv
* 007      G. Pellerin (Feb 96) Revised calculation for H (stable)
* 008      G. Pellerin (Feb 96) Remove corrective terms to CTU
* 009      Y. Delage and B. Bilodeau (Jul 97) - Cleanup
* 010      Y. Delage (Feb 98) - Addition of HMIN
* 011      Y. Delage (Sept 00) - Set top of surface layer at ZU +Z0
*                              - Output UE instead of UE**2
*                              - Initialise ILMO and H
*                              - Change iteration scheme for stable case
*                              - Introduce log-linear profile for near-
*                                 neutral stable cases
*                              - set VAMIN inside flxsurf
* 012      Y. Delage (Oct 00) - Input of wind and temperatur/humidity
*                                at different levels
*
*Object
*          to calculate surface layer transfer coefficients and fluxes
*          FLXSURFZ is a variant of FLXSURF3 that permits to input
*          wind and temperature (humidity) at different levels
*
*Arguments
*
*          - Output -
* CMU      transfer coefficient of momentum times UE
* CTU      transfer coefficient of temperature times UE
* RIB      bulk Richardson number
* FTEMP    temperature flux
* FVAP     vapor flux
* ILMO     (1/length of Monin-Obukov)
* UE       friction velocity
* H        height of the boundary layer
* FM       momentum stability function
* FH       heat stability function
* LZZ0     log ((zu+z0)/z0)
* LZZ0T    log ((zt+z0)/z0t)
*          - Input -
* FCOR     Coriolis factor
* ZU       height of wind input
* ZT       height of temperature and humidity input
* TA       temperature  at z=zt
* QA       specific humidity at z=zt
* VA       wind speed at z=zu
* TG       surface temperature
* QG       specific humidity at the surface
* Z0       roughness length for momentum      flux calculations
* Z0T      roughness length for heat/moisture flux calculations
* N        horizontal dimension
* ITER     key allowing to bypass the calculations
*
*Notes
*          SEE DELAGE AND GIRARD BLM 58 (19-31)
*                "       BLM 82 (23-48)
*
*     DIVERSES CONSTANTES PHYSIQUES
*
c
      REAL CPD, CPV, RGASD, RGASV, TRPL, TCDK, RAUW, EPS1
      REAL DELTA, CAPPA, TGL, CONSOL, GRAV, RAYT, STEFAN, PI
      REAL OMEGA, EPS2
      REAL KNAMS, STLO, KARMAN, RIC, CHLC, CHLF, DG2RAD
*
      COMMON/CTSPHY/  CPD, CPV, RGASD, RGASV, TRPL, TCDK, RAUW,
     $                EPS1, EPS2, DELTA, CAPPA, TGL, CONSOL,
     $                GRAV, RAYT, STEFAN, PI, OMEGA,
     $                KNAMS, STLO, KARMAN, RIC, CHLC, CHLF, DG2RAD
      REAL AS,ASX,CI,BS,BETA,FACTN,HMIN
c
      COMMON /SURFCON/AS,ASX,CI,BS,BETA,FACTN,HMIN
      REAL TPOTA(N)
      INTEGER J
      INTEGER IT,ITMAX
      REAL HMAX,CORMIN,EPSLN
      REAL RAC3,CM,CT,ZP
      REAL F,GG,DG,X,X0,Y,Y0
      REAL UNSH,HE,HS
      REAL DTHV,TVA,TVS
      REAL HL,VAMIN,U
      REAL CS
      REAL ZB,DD,ILMOX
      REAL DF,ZZ
      SAVE HMAX,CORMIN,EPSLN
      SAVE ITMAX
      SAVE VAMIN
*
      DATA CORMIN, HMAX /0.7E-4 ,  1500.0/
      DATA ITMAX / 3  /
      DATA VAMIN /0.1/
      DATA EPSLN / 1.0e-05 /
*
      real a,b,c,d,psi,z
      real unsl,hi
************************************************************************
**  fonctions de couche de surface pour le cas stable                 **
************************************************************************
*
      d  (unsl) = 4*AS*BETA*unsl
      c  (hi)   = d(unsl)*hi - hi**2
      b  (hi)   = d(unsl) - 2*hi
      a  (z,hi) = sqrt(1 + b(hi)*z - c(hi)*z**2)
      psi(z,hi) = 0.5 * (a(z,hi)-z*hi-log(1+b(hi)*z*0.5+a(z,hi))-
     +            b(hi)/(2*sqrt(c(hi)))*asin((b(hi)-2*c(hi)*z)/d(unsl)))
*
*   Limites de validite: unsl >= 0 (cas stable ou neutre)
*                        c > 0 (hi < d)
*                        z*hi < 1
*   Ces 2 dernieres conditions imposees a l'aide du facteur 'factn'
*
*   Reference :  Y. Delage, BLM, 82 (p23-48) (Eq.33-37)
************************************************************************
      DF(ZZ)=(1-ZZ*UNSH)*sqrt(1+d(ilmo(j))*ZZ/(1-ZZ*UNSH))
*
*     CONVERTIR TA A TEMPERATURE POTENTIELLE TPOTA
      DO J=1,N
         TPOTA(J) = TA(J) + (GRAV/CPD)*ZT(J)
      END DO
*
      RAC3=sqrt(3.0)
      CS=AS*2.5
*
      DO 1 J=1,N
      IF(ITER(J).EQ.1) THEN
*
*  CALCULATE THE RICHARDSON NUMBER
        ZP=ZU(J)**2/(ZT(J)+Z0(J)-Z0T(J))
        u=max(vamin,REAL(va(j)))
        tva=(1+DELTA*QA(J))*TPOTA(J)
        tvs=(1+DELTA*QG(J))*TG(J)
        dthv=tva-tvs
        RIB(J)=grav/(tvs+0.5*dthv)*ZP*dthv/u**2
        if (rib(j).ge.0.0) rib(j) = max(rib(j), EPSLN)
        if (rib(j).lt.0.0) rib(j) = min(rib(j),-EPSLN)
*
*  FIRST APPROXIMATION TO ILMO
        LZZ0(J)=LOG(1+ZU(J)/Z0(J))
        LZZ0T(J)=LOG((ZT(J)+Z0(J))/Z0T(J))
        IF(RIB(J).GT.0.)  THEN
           FM(J)=LZZ0(J)+CS*RIB(J)/max(REAL(2*z0(j)),1.0)
           FH(J)=BETA*(LZZ0T(J)+CS*RIB(J))/
     1           max(REAL(sqrt(z0(j)*z0t(j))),1.0)
           ILMO(J)=RIB(J)*FM(J)**2/(ZP*FH(J))
           F=MAX(REAL(ABS(FCOR(J))),CORMIN)
           H(J)=BS*sqrt(karman*u/(ILMO(J)*F*fm(j)))
        ELSE
           FM(J)=LZZ0(J)-min(0.7+log(1-rib(j)),LZZ0(J)-1)
           FH(J)=BETA*(LZZ0T(J)-min(0.7+log(1-rib(j)),LZZ0T(J)-1))
           ILMO(J)=RIB(J)*FM(J)**2/(ZP*FH(J))
        ENDIF
      ENDIF
    1 CONTINUE
* - - - - - - - - -  BEGINNING OF ITERATION LOOP - - - - - - - - - - -
      DO 35 IT=1,ITMAX
      DO 35 J=1,N
      IF(ITER(J).EQ.1) THEN
        u=max(vamin,REAL(va(j)))
        ZP=ZU(J)**2/(ZT(J)+Z0(J)-Z0T(J))
        IF(RIB(J).GT.0.)  THEN
*----------------------------------------------------------------------
*  STABLE CASE
       	ILMO(J)=max(EPSLN,ILMO(J))
        hl=(ZU(J)+10*Z0(J))*FACTN
        F=MAX(REAL(ABS(FCOR(J))),CORMIN)
        hs=BS*sqrt(karman*u/(ILMO(J)*F*fm(j)))
        H(J)=MAX(HMIN,hs,hl,factn/d(ILMO(J)))
        UNSH=1/H(J)
        unsl=ILMO(J)
        fm(J)=LZZ0(J)+psi(REAL(ZU(J)+Z0(J)),unsh)-psi(REAL(Z0(J)),unsh)
        fh(J)=BETA*(LZZ0T(J)+psi(REAL(ZT(J)+Z0(J)),unsh)
     *       -psi(REAL(Z0T(J)) ,unsh))
        DG=(-ZP)*FH(J)/FM(J)**2*(1+beta*(DF(REAL(ZT(J)+Z0(J)))
     *       -DF(REAL(Z0T(J))))/
     1          (2*FH(J))-(DF(REAL(ZU(J)+Z0(J)))-DF(REAL(Z0(J))))/FM(J))

*----------------------------------------------------------------------
*  UNSTABLE CASE
      ELSE
        ILMO(J)=MIN(0.,ILMO(J))
         X=(1-CI*(ZU(J)+Z0(J))*BETA*ILMO(J))**(1./6)
         X0=(1-CI*Z0(J)*BETA*ILMO(J))**(1./6.)
         FM(J)=LZZ0(J)+LOG((X0+1)**2*SQRT(X0**2-X0+1)*(X0**2+X0+1)**1.5
     %               /((X+1)**2*SQRT(X**2-X+1)*(X**2+X+1)**1.5))
     %              +RAC3*ATAN(RAC3*((X**2-1)*X0-(X0**2-1)*X)/
     %              ((X0**2-1)*(X**2-1)+3*X*X0))
         Y=(1-CI*(ZT(J)+Z0(J))*BETA*ILMO(J))**(1./3)
         Y0=(1-CI*Z0T(J)*BETA*ILMO(J))**(1./3)
         FH(J)=BETA*(LZZ0T(J)+1.5*LOG((Y0**2+Y0+1)/(Y**2+Y+1))+RAC3*
     %        ATAN(RAC3*2*(Y-Y0)/((2*Y0+1)*(2*Y+1)+3)))
         DG=(-ZP)*FH(J)/FM(J)**2*(1+beta/FH(J)*(1/Y-1/Y0)-2/FM(J)*
     %                (1/X-1/X0))
      ENDIF
*----------------------------------------------------------------------
      IF(IT.LT.ITMAX) THEN
             GG=RIB(J)-FH(J)/FM(J)**2*ZP*ILMO(J)
             ILMO(J)=ILMO(J)-GG/DG
      ENDIF
      ENDIF
   35 CONTINUE
* - - - - - -  - - - END OF ITERATION LOOP - - - - - - - - - - - - - -
      DO 80 J=1,N
      IF(ITER(J).EQ.1) THEN
        u=max(vamin,REAL(va(j)))
*----------------------------------------------------------------------
*  CALCULATE ILMO AND STABILITY FUNCTIONS FROM LOG-LINEAR PROFILE
        zb=zu(j)/(zt(j)+z0(j)-z0t(j))
        dd=(beta*lzz0t(j)*zb)**2-4*rib(j)*asx*lzz0(j)*beta*
     1        (lzz0t(j)*zb-lzz0(j))
        if(rib(j).gt.0..and.rib(j).lt.beta/asx.and.dd.ge.0.) then
         ilmox=(beta*lzz0t(j)*zb-2*rib(j)*asx*lzz0(j)-sqrt(dd))
     1          /(2*zu(j)*(asx**2*rib(j)-beta*asx))
         if(ilmox.lt.ilmo(j)) then
            ilmo(j)=ilmox
            fm(j)=lzz0(j)+asx*zu(j)*ilmox
            fh(j)=beta*(lzz0t(j)+asx*(zt(j)+z0(j)-z0t(j))*ilmox)
         endif
        endif
*----------------------------------------------------------------------
        CM=karman/FM(J)
        CT=karman/FH(J)
        UE(J)=u*CM
        CMU(J)=CM*UE(J)
        CTU(J)=CT*UE(J)
        if (rib(j).gt.0.0) then
*             cas stable
              H(J)=MIN(H(J),hmax)
        else
*             cas instable
              F=MAX(ABS(REAL(FCOR(J))),CORMIN)
              he=max(HMIN,0.3*UE(J)/F)
              H(J)=MIN(he,hmax)
        endif
        FTEMP(J)=(-CTU(J))*(TPOTA(J)-TG(J))
        FVAP(J)=(-CTU(J))*(QA(J)-QG(J))
      ENDIF
   80 CONTINUE
      RETURN
      END

max_stars_repo_name: OSADP/Pikalert-Vehicle-Data-Translator-

max_issues_repo_name: OSADP/Pikalert-Vehicle-Data-Translator-

max_forks_repo_name: OSADP/Pikalert-Vehicle-Data-Translator-


text:
! Calculate SUITS coefficients
subroutine SUITS(na,litab,lidf,tts,tto,cts,cto,psi,ctscto, &
        ks,ko,sob,sof,sdb,sdf,dob,dof,ddb,ddf)

    implicit none
    integer, intent(in) :: na
    REAL*8, INTENT(in) :: litab(na),lidf(na),tts,tto,cts,cto,psi,ctscto
    REAL*8, INTENT(out) :: ks,ko,sob,sof,sdb,sdf,dob,dof,ddb,ddf

    real*8 :: pi,rd
    common /angle/ pi,rd
    REAL*8 :: bf,ttl,ctl,chi_s,chi_o,frho,ftau,ksli,koli,sobli,sofli,bfli
    integer :: i

    !	Calculate geometric factors associated with extinction and scattering
    !	Initialise sums
    ks  = 0
    ko  = 0
    bf  = 0
    sob = 0
    sof = 0

    !	Weighted sums over LIDF
    DO i = 1,na
        ttl = litab(i)      ! leaf inclination discrete values
        ctl = COS(rd*ttl)
        !	SAIL volume scattering phase function gives interception and portions to be
        !	multiplied by rho and tau
        CALL volscatt(tts,tto,psi,ttl,chi_s,chi_o,frho,ftau)

        !********************************************************************************
        !*                   SUITS SYSTEM COEFFICIENTS
        !*
        !*	ks  : Extinction coefficient for direct solar flux
        !*	ko  : Extinction coefficient for direct observed flux
        !*	att : Attenuation coefficient for diffuse flux
        !*	sigb : Backscattering coefficient of the diffuse downward flux
        !*	sigf : Forwardscattering coefficient of the diffuse upward flux
        !*	sf  : Scattering coefficient of the direct solar flux for downward diffuse flux
        !*	sb  : Scattering coefficient of the direct solar flux for upward diffuse flux
        !*	vf   : Scattering coefficient of upward diffuse flux in the observed direction
        !*	vb   : Scattering coefficient of downward diffuse flux in the observed direction
        !*	w   : Bidirectional scattering coefficient
        !********************************************************************************

        !	Extinction coefficients
        ksli = chi_s/cts
        koli = chi_o/cto

        !	Area scattering coefficient fractions
        sobli = frho*pi/ctscto
        sofli = ftau*pi/ctscto
        bfli  = ctl*ctl
        ks    = ks+ksli*lidf(i)
        ko    = ko+koli*lidf(i)
        bf    = bf+bfli*lidf(i)
        sob   = sob+sobli*lidf(i)
        sof   = sof+sofli*lidf(i)

    ENDDO

    !	Geometric factors to be used later with rho and tau
    sdb = 0.5*(ks+bf)
    sdf = 0.5*(ks-bf)
    dob = 0.5*(ko+bf)
    dof = 0.5*(ko-bf)
    ddb = 0.5*(1.+bf)
    ddf = 0.5*(1.-bf)

    return
end subroutine SUITS

! *****************

!! Correct rho and tau using SUITS factors
subroutine RTgeom(rho,tau,ddb,ddf,sdb,sdf,dob,dof,sob,sof, &
        sigb,att,m,sb,sf,vb,vf,w)

    use mod_dataSpec_wavelength
    implicit none
    real*8,intent(in),dimension(nw) :: rho,tau
    real*8,intent(in) :: ddb,ddf,sdb,sdf,dob,dof,sob,sof
    real*8,intent(out),dimension(nw) :: sigb,att,m,sb,sf,vb,vf,w

    real*8,dimension(nw) :: m2,sigf

    sigb = ddb*rho+ddf*tau
    sigf = ddf*rho+ddb*tau
    att  = 1.-sigf
    m2   = (att+sigb)*(att-sigb)
    WHERE (m2 < 0)
        m2 = 0
    ENDWHERE
    m  = SQRT(m2)
    sb = sdb*rho+sdf*tau
    sf = sdf*rho+sdb*tau
    vb = dob*rho+dof*tau
    vf = dof*rho+dob*tau
    w  = sob*rho+sof*tau

    return
end subroutine RTgeom

max_stars_repo_licenses: NCSA Unlicense

max_issues_repo_licenses: NCSA Unlicense

max_forks_repo_licenses: NCSA Unlicense


text:
PROGRAM me1xb

! Code converted using TO_F90 by Alan Miller
! Date: 2004-06-28  Time: 12:58:10

!       =========================================================
!       Purpose: This program computes the exponential integral
!                E1(x) using subroutine E1XB
!       Input :  x  --- Argument of E1(x)  ( x > 0 )
!       Output:  E1 --- E1(x)
!       Example:
!                  x          E1(x)
!                -------------------------
!                 0.0     .1000000000+301
!                 1.0     .2193839344E+00
!                 2.0     .4890051071E-01
!                 3.0     .1304838109E-01
!                 4.0     .3779352410E-02
!                 5.0     .1148295591E-02
!       =========================================================

DOUBLE PRECISION :: e1,x
WRITE(*,*)'Please enter x '
!        READ(*,*) X
x=5.0
WRITE(*,*)'   x          E1(x)'
WRITE(*,*)' -------------------------'
CALL e1xb(x,e1)
WRITE(*,10)x,e1
10      FORMAT(1X,f5.1,e20.10)
END PROGRAM me1xb

SUBROUTINE e1xb(x,e1)

!       ============================================
!       Purpose: Compute exponential integral E1(x)
!       Input :  x  --- Argument of E1(x)
!       Output:  E1 --- E1(x)  ( x > 0 )
!       ============================================

DOUBLE PRECISION, INTENT(IN)             :: x
DOUBLE PRECISION, INTENT(OUT)            :: e1
IMPLICIT DOUBLE PRECISION (a-h,o-z)
IF (x == 0.0) THEN
  e1=1.0D+300
ELSE IF (x <= 1.0) THEN
  e1=1.0D0
  r=1.0D0
  DO  k=1,25
    r=-r*k*x/(k+1.0D0)**2
    e1=e1+r
    IF (DABS(r) <= DABS(e1)*1.0D-15) EXIT
  END DO
  15         ga=0.5772156649015328D0
  e1=-ga-DLOG(x)+x*e1
ELSE
  m=20+INT(80.0/x)
  t0=0.0D0
  DO  k=m,1,-1
    t0=k/(1.0D0+k/(x+t0))
  END DO
  t=1.0D0/(x+t0)
  e1=DEXP(-x)*t
END IF
RETURN
END SUBROUTINE e1xb

max_stars_repo_name: vasaantk/bin

max_issues_repo_name: vasaantk/bin

max_forks_repo_name: vasaantk/bin