text:
!> \file covariance.f90
!! @todo rewrite each program/module that uses the old module \a  covariance to use the \a evd_covariance module. The \a evd_covariance module works together with \a kronecker modules
!!
!! \brief General module for covariance matrix
!!
!! @author Innocent Souopgui
!!
! !! \note In recent update, the old covariance module was separated into Kronecker module and the current covariance module. User need to reference the kronecker module for compilation
! !!
! !! \details This module implements a set of subroutines for the covariance matrix for one and two dimensional problems
! !! Using the positive (semi) definiteness of the covariance matrices, the covariance matrix B is diagonalizable
! !! and its singular value decomposition is also its diagonalization in a basis of eigenvectors.
! !! We consider the case where B is approximated using only the leading eigenvectors, the full case being similar.
! !! The covariance matrix is supposed to be given under a factorized forms
! !!   \f[ B \approx UDU^T \f] where the colums of U are the leading eigenvectors
! !!   of B and D is the diagonal matrix of the leading eigenvalues of B.
! !! For the 2D, we consider the covariance function that factorises as a product of covariance functions
! !! in each spatial dimension; so that the covariance matrix B is the kronecker product of the covariance
! !! matrices in the single dimensions \f[ B = Kron(Bx,By) \f]. Bx and By are given under their factorized forms
! !! \f[ B Kron(Ux,Uy)*Kron(Dx,Dy)*Kron(Ux,Uy)^T \f], where Dx is a diagonal matrix with leading eigenvalues
! !! of Bx on the diagonal and Dy is a diagonal matrix with eigenvalues of By on the diagonal
! !! Bx = Ux*Dx*Ux^T,  By = Uy*Dy*Uy^T.
! !!The 2D covariance matrix applies to a vector that correspond to a discretization of
! !!  a 2d function (function of two variables) on the same grid as the control
! !!  variable. The chronecker product factorization depends on the ordering of
! !!  the dimensions in the discretization as well as the vectorization of the
! !!  2D discretized function.
! !!Suppose that the covariance function factorizes such that the covariances on
! !!  the grid are combinations of the covariances in each direction,
! !!  cov(p1, p2) = cov_x(x1,x2)cov_y(y1,y2) where p1 is the point of coordinates
! !!  (x1,y1) and p2 (x2,y2). In this case B*f = Kron(B_x,B_y)*f where B_x and
! !!  B_y are respectively the covariances in the direction x, resp. y.
! !! According to the ordering of the discretization points and the vectorization
! !!   order, that formula may change:
! !! Let F be the matrix of a discretized function, and f its vectorized representation.
! !!if the first dimension of F corresponds to the x direction (x-leading) then,
! !!  g=B*f = Kron(B_x,B_y)*f, if f is a row major vectorization (C-prefered ordering);
! !!  and g=B*f = Kron(B_y,B_x)*f, if f is a column major vectorization
! !!  (Fortran-prefered ordering).
! !!if the first dimension of F corresponds to the y direction (y-leading) then,
! !! g=B*f = Kron(B_y,B_x)*f with a row major vectorization (C-prefered ordering)
! !! of F; and g=B*f = Kron(B_x,B_y)*f, with a column major vectorization
! !! (Fortran-prefered ordering) of F.
! !! All the 2D routines in this module supposed either:
! !! x-leading and and row major vectorization or
! !! y-leading and the column major vectorization
! !!
! !<
! MODULE covariance
!   USE general_constant
!   USE general_tools, only: readInfo, readMatrix
!   USE debug_tools
!   use kronecker
! IMPLICIT NONE
!
!   private
!   public svd_cov_matrix, load_svd_covariance, apply_svd_cov&
!     , apply_svd_cov_inv, apply_svd_square_cov, get_svd_cov_size&
!     , apply_svd_square_covb, apply_svd_square_cov_inv, get_svd_cov_nsv
!   !, cov_mv, cov_inv_mv, cov_square_mv, cov_square_mvb, cov_square_inv_mv&
!
!
!   !> \brief default number of factors in the Kronecker factorization of the covariance matrix, used for non initialized covariance matrix
!   !<
!   INTEGER, PARAMETER :: default_nFact=-1
!
! !> Derived type for singular values decomposed covariance matrices
! !! It is assumed that the covariance matrix is provided as its Singular Value Decomposition: B = UDUt. Furthermore, it can be (Kronecker) factorized along the physical dimensions of the problem: B = Bx kron By kron Bz
! type svd_cov_matrix
!
!   !> \brief number of terms in the Kronecker factorization of the covariance matrix.
!   !!Conceptually, it corresponds to the number of physical dimensions of the problem if the covariance is factorized.
!   !!Logically, it can be any Kronecker factorization.
!   !<
!   INTEGER :: im_nFac = default_nFact
!
!   !> \brief matrix(ces) of singular vector of the covariance matrix
!   !! If the covariance matrix B is not factorized, only Ux is initialized
!   !! IF B is factorized as B = Bx kron By (two terms), Ux and Uy are provided
!   !! A maximum of 3 terms in the factorization is supported. Only 1 and 2 terms are implemented for now.
!   !<
!   REAL(cp), ALLOCATABLE, DIMENSION(:,:):: Ux, Uy, Uz
!
!   !> vector(s) of singular values of the SVD of B
!   REAL(cp), ALLOCATABLE, DIMENSION(:) :: Dx, Dy, Dz
!
! end type svd_cov_matrix
!
! !private procedures interface
!   !> Covariance matrix vector product
!   INTERFACE cov_mv
!     MODULE PROCEDURE cov_mv_default
!     MODULE PROCEDURE cov_mv_2Fac
!     !MODULE PROCEDURE cov_mv_3Fac
!   END INTERFACE cov_mv
!
!   !> Inverse covariance matrix vector product
!   INTERFACE cov_inv_mv
!     MODULE PROCEDURE cov_inv_mv_default
!     MODULE PROCEDURE cov_inv_mv_2Fac
!     !MODULE PROCEDURE cov_inv_mv_3Fac
!   END INTERFACE cov_inv_mv
!
!   !> square root covariance matrix vector product
!   INTERFACE cov_square_mv
!     MODULE PROCEDURE cov_square_mv_default
!     !MODULE PROCEDURE cov_square_mv_1D
!     MODULE PROCEDURE cov_square_mv_2Fac
!   END INTERFACE cov_square_mv
!
!   !> square root covariance matrix vector product (adjoint)
!   INTERFACE cov_square_mvb
!     MODULE PROCEDURE cov_square_mvb_default
!     MODULE PROCEDURE cov_square_mvb_2Fac
!     !MODULE PROCEDURE cov_square_mv_3Fac
!   END INTERFACE cov_square_mvb
!
!   !> Inverse square root covariance matrix vector product
!   INTERFACE cov_square_inv_mv
!     MODULE PROCEDURE cov_square_inv_mv_default
!     MODULE PROCEDURE cov_square_inv_mv_2Fac
!     !MODULE PROCEDURE cov_square_inv_mv_3Fac
!   END INTERFACE cov_square_inv_mv
!
!   !> Checks the parameters for the application of the covariance matrix or its inverse
!   INTERFACE cov_mv_check
!     MODULE PROCEDURE cov_mv_check_default
!     MODULE PROCEDURE cov_mv_check_2Fac
!   END INTERFACE cov_mv_check
!
!   !> Check the parameters for the application of the square root of the covariance matrix or its inverse
!   INTERFACE cov_square_mv_check
!     MODULE PROCEDURE cov_square_mv_check_default
!     MODULE PROCEDURE cov_square_mv_check_2Fac
!   END INTERFACE cov_square_mv_check
! !end of private procedures interface
!
! CONTAINS
!
! !public procedures
!
!   !> \brief returns the number of singular values used to define the approximation of the covariance matrix; that is also the size of the changed state variable
!   !! \param[in] C covariance matrix derived type
!   !! \param[in] [nctl] size of the control vector
!   !! If a call to load_covariance has not been done before: the function returns nctl if it is provided, otherwise, the program is stopped. If a call to load_covariance have been done before, the associated value is returned
!   !<
!   FUNCTION get_svd_cov_nsv(C, nctl) RESULT(nsv)
!     type(svd_cov_matrix), intent(in) :: C
!     INTEGER, OPTIONAL, INTENT(IN) :: nctl
!     !local variables
!     INTEGER :: nsv
!
!     SELECT CASE(C%im_nFac)
!       CASE (1)
!         nsv = SIZE(C%Dx)
!       CASE (2)
!         nsv = SIZE(C%Dx)*SIZE(C%Dy)
!       CASE (default_nFact)!default value
!         IF( PRESENT(nctl) )THEN
!           nsv = nctl
!         ELSE
!           CALL stop_program( 'In svd_cov_nsv, covariance not loaded and nctl not provided' )
!         END IF
!       CASE DEFAULT
!         CALL stop_program( 'In svd_cov_nsv, covariance operations are define only for 1D and 2D' )
!     END SELECT
!   END FUNCTION get_svd_cov_nsv
!
!   !> \brief Reads a covariance matrix stored under the form of an SVD decomposition
!   !! \param[in] UFName file name of the matrix of singular vectors
!   !! \param[in] DFName file name of the matrix of singular values
!   !! \param[inout] U matrix of singular vectors
!   !! \param[inout] D vector of singular values
!   !<
!   SUBROUTINE readSVD( UFName, DFName, U, D )
!     CHARACTER(len=*), INTENT(IN)  :: UFName, DFName
!     REAL(cp), DIMENSION(:,:), INTENT(INOUT) :: U
!     REAL(cp), DIMENSION(:), INTENT(INOUT) :: D
!     !local variables
!     REAL(cp), DIMENSION(:,:), ALLOCATABLE :: T
!     INTEGER :: il_n, il_k, il_nRow, il_nCol, ibj
!
!     CALL readSVDSize( UFName, DFName, il_n, il_k, il_nRow, il_nCol )
!     ALLOCATE( T(il_nRow, il_nCol) )
!
!     CALL readMatrix( UFName, U )
!     CALL readMatrix( DFName, T )
!
!     IF(il_nRow == 1) THEN
!       D = T(1,:)
!     ELSE IF(il_nCol == 1) THEN
!       D = T(:,1)
!     ELSE
!       DO ibj = 1, il_nCol
!         D(ibj) = T( ibj, ibj )
!       END DO
!     END IF
!   END SUBROUTINE readSVD
!
!   !> \brief read the size of the SVD decomposition of B: size of each singular vector, and number of singular vector
!   !! \param[in] UFName file name of the matrix of singular vectors
!   !! \param[in] DFName file name of the matrix of singular values
!   !! \param[inout] n size of each singular vector (number of rows/colums of B)
!   !! \param[inout] k number of singular values (or vector)
!   !! \param[inout] nrD (optional) number of rows of D
!   !! \param[inout] ncD (optional) number of colums of D
!   !! if the sizes in UFName and DFName do not match, the program terminates
!   !<
!   SUBROUTINE readSVDSize(UFName, DFName, n, k, nrD, ncD)
!     CHARACTER(len=*), INTENT(IN)  :: UFName, DFName
!     INTEGER, INTENT(OUT) :: n, k
!     INTEGER, OPTIONAL, INTENT(OUT) :: nrD, ncD
!     INTEGER :: il_nrD, il_ncD
!
!     CALL readInfo( UFName, n, k )
!     CALL readInfo( DFName, il_nrD, il_ncD )
!     CALL SVD_check( n, k, il_nrD, il_ncD )
!     IF( PRESENT(nrD) ) nrD = il_nrD
!     IF( PRESENT(ncD) ) ncD = il_ncD
!   END SUBROUTINE readSVDSize
!
!   !> \brief Loads the covariance matrix represented by its singular value decomposition
!   !! \param[in, out] C covariance matrix derived type
!   !! \param[in] dir directory containing the data files
!   !!
!   !! \details It is assumed that the covariance matrix B is provided as its
!   !!  Singular Value Decomposition: B = UDU^t. Furthermore, it can be
!   !!  (Kronecker) factorized (for example along the physical dimensions of
!   !!  the problem) in which case it is given under the form:
!   !! - B = Bx kron By = (UxDxUx^t) kron (UyDyUy^t) or
!   !! - B = Bx kron By kron Bz = (UxDxUx^t) kron (UyDyUy^t) kron (UzDzUz^t)
!   !! Conventions for the file names. Data files are named after the matrices names with the extension 'dat'
!   !! - 'U.dat' and 'D.dat' for the non factorized case
!   !! - 'Ux.dat',  'Uy.dat', 'Uz.dat' and 'Dx.dat', 'Dy.dat', 'Dz.dat' for the factorized case
!   !! the existence of files is checked to determine the appropriate configuration.
!   !! ('U.dat' and 'D.dat') or ('Ux.dat' and 'Dx.dat') should be present, otherwise the program terminates
!   !!
!   !<
!   SUBROUTINE load_svd_covariance(C, dir)
!     type(svd_cov_matrix), intent(in out) :: C
!     CHARACTER(len=*), INTENT(IN)  :: dir
!     !local variables
!     INTEGER :: il_N, il_K
!     CHARACTER(len=ip_fnl) :: UFName, DFName
!
!     CALL debug('Entering load_covariance ', tag=dALLWAYS )
!     !Checking files
!     IF( ACCESS(TRIM(dir)//'/U.dat', MODE='r')==0 )THEN
!       UFName = TRIM(dir)//'/U.dat'
!       DFName = TRIM(dir)//'/D.dat'
!     ELSE
!       UFName = TRIM(dir)//'/Ux.dat'
!       DFName = TRIM(dir)//'/Dx.dat'
!     END IF
!
!     CALL debug(100, ' In load_svd_covariance ', tag=dALLWAYS )
!     CALL readSVDSize(UFName, DFName, il_N, il_K)
!     ALLOCATE( C%Ux(il_N, il_K), C%Dx(il_K) )
!     CALL readSVD( UFName, DFName, C%Ux, C%Dx )
!     CALL debug(C%Ux, 'Ux = ')
!     CALL debug(C%Dx, 'Dx = ')
!     C%im_nFac = 1
!
!     CALL debug(200, ' In load_svd_covariance ', tag=dALLWAYS )
!     UFName = TRIM(dir)//'/Uy.dat'
!     DFName = TRIM(dir)//'/Dy.dat'
!     IF( ACCESS(uFName, MODE='r')==0 )THEN
!       CALL readSVDSize(UFName, DFName, il_N, il_K)
!       ALLOCATE( C%Uy(il_N, il_K), C%Dy(il_K) )
!       CALL readSVD( UFName, DFName, C%Uy, C%Dy )
!       CALL debug(C%Uy, 'Uy = ')
!       CALL debug(C%Dy, 'Dy = ')
!       C%im_nFac = C%im_nFac + 1
!     END IF
!
!     CALL debug('Exiting load_svd_covariance ', tag=dALLWAYS )
!   END SUBROUTINE load_svd_covariance
!
!   !> Gets the sizes of the factorized arrays that make up the EV decomposed covariance matrix
!   !! \param[in] C covariance matrix derived type
!   !! \param[out] Nx number of components [in the x direction]
!   !! \param[out] Kx number of singular vectors [in the x direction]
!   !! \param[out] Ny number  of components in the y direction
!   !! \param[out] Ky number of singular vectors in the y direction
!   !!
!   !! \details Negatives value means that there is no component in the given direction.
!   !<
!   subroutine get_svd_cov_size(C, Nx, Kx, Ny, Ky)
!     type(svd_cov_matrix), intent(in out) :: C
!     INTEGER, INTENT(OUT), OPTIONAL :: Nx, Kx, Ny, Ky
!
!     if(C%im_nFac>=1)then
!       if( present(Nx) )Nx = size(C%Ux,1)
!       if( present(Kx) )Kx = size(C%Ux,2)
!     else
!       if( present(Nx) )Nx = -1
!       if( present(Kx) )Kx = -1
!     end if
!
!     if(C%im_nFac>=2)then
!       if( present(Ny) )Ny = size(C%Uy,1)
!       if( present(Ky) )Ky = size(C%Uy,2)
!     else
!       if( present(Ny) )Ny = -1
!       if( present(Ky) )Ky = -1
!     end if
!   end subroutine get_svd_cov_size
!
!   !> \brief applies the covariance matrix to the vector f, g = Bf
!   !! \param[in] C covariance matrix derived type
!   !! \param[in] f input vector
!   !! \param[in] g output vector
!   !! \param[in] [die] specify the behavior of the program if the covariance matrix is not loaded. if TRUE, the program stops, otherwise g =f. The user should make sure that f and g have the same size.
!   !<
!   SUBROUTINE apply_svd_cov(C, f, g, die)
!     type(svd_cov_matrix), intent(in) :: C
!     REAL(cp), DIMENSION(:), INTENT(IN)  :: f
!     REAL(cp), DIMENSION(:), INTENT(OUT) :: g
!     LOGICAL, OPTIONAL, INTENT(IN) :: die
!     !local variables
!     LOGICAL :: ll_die
!
!     IF( PRESENT(die) )THEN
!       ll_die = die
!     ELSE
!       ll_die = .TRUE.
!     END IF
!
!     SELECT CASE(C%im_nFac)
!       CASE (1)
!         CALL stop_program('apply_svd_cov is not yed defined for 1D')
!       CASE (2)
!         CALL cov_mv(C%Ux, C%Dx, C%Uy, C%Dy, f, g)
!       CASE (default_nFact)!default value
!         IF(ll_die)THEN
!           CALL stop_program('In apply_svd_cov covariance matrix is not loaded')
!         ELSE
!           g=f
!         END IF
!     END SELECT
!   END SUBROUTINE apply_svd_cov
!
!   !> \brief applies the inverse of covariance matrix, g = B^-1 f
!   !! \param[in] C covariance matrix derived type
!   !! \param[in] f input vector
!   !! \param[in] g output vector
!   !! \param[in] [die] specify the behaviour of the program if the covariance matrix is not loaded. if TRUE, the program stops, otherwise g =f. The user should make sure that f and g have the same size.
!   !<
!   SUBROUTINE apply_svd_cov_inv(C, f, g, die)
!     type(svd_cov_matrix), intent(in) :: C
!     REAL(cp), DIMENSION(:), INTENT(IN)  :: f
!     REAL(cp), DIMENSION(:), INTENT(OUT) :: g
!     LOGICAL, OPTIONAL, INTENT(IN) :: die
!     !local variables
!     LOGICAL :: ll_die
!
!     IF( PRESENT(die) )THEN
!       ll_die = die
!     ELSE
!       ll_die = .TRUE.
!     END IF
!
!     SELECT CASE(C%im_nFac)
!       CASE (1)
!         CALL stop_program('apply_svd_cov_inv is not yed defined for 1D')
!       CASE (2)
!         CALL cov_inv_mv(C%Ux, C%Dx, C%Uy, C%Dy, f, g)
!       CASE (default_nFact)!default value
!         IF(ll_die)THEN
!           CALL stop_program('In apply_svd_cov_inv covariance matrix is not loaded')
!         ELSE
!           g=f
!         END IF
!     END SELECT
!   END SUBROUTINE apply_svd_cov_inv
!
!   !> \brief apply the square root of the covariance matrix to the vector f and result g
!   !! \param[in] C covariance matrix derived type
!   !! \param[in] f input vector
!   !! \param[in] g output vector
!   !! \param[in] [die] specify the behaviour of the program if the covariance matrix is not loaded. if TRUE, the program stops, otherwise g =f, the user should make sure that f and g have the same size.
!   !<
!   SUBROUTINE apply_svd_square_cov(C, f, g, die)
!     type(svd_cov_matrix), intent(in) :: C
!     REAL(cp), DIMENSION(:), INTENT(IN)  :: f
!     REAL(cp), DIMENSION(:), INTENT(OUT) :: g
!     LOGICAL, OPTIONAL, INTENT(IN) :: die
!     !local variables
!     LOGICAL :: ll_die
!
!     IF( PRESENT(die) )THEN
!       ll_die = die
!     ELSE
!       ll_die = .TRUE.
!     END IF
!
!     SELECT CASE(C%im_nFac)
!       CASE (1)
!         CALL stop_program('apply_svd_square_cov is not yed defined for 1D')
!       CASE (2)
!         CALL cov_square_mv(C%Ux, C%Dx, C%Uy, C%Dy, f, g)
!       CASE (default_nFact)!default value
!         IF(ll_die)THEN
!           CALL stop_program('In apply_svd_square_cov, covariance matrix is not loaded')
!         ELSE
!           g=f
!         END IF
!     END SELECT
!   END SUBROUTINE apply_svd_square_cov
!
!   !> \brief apply the adjoint of the square root of the covariance matrix to the vector f and result g
!   !! \param[in] C covariance matrix derived type
!   !! \param[in] f input vector
!   !! \param[in] g output vector
!   !! \param[in] [die] specify the behaviour of the program if the covariance matrix is not loaded. if TRUE, the program stops, otherwise g =f. The user should make sure that f and g have the same size.
!   !<
!   SUBROUTINE apply_svd_square_covb(C, fb, gb, die)
!     type(svd_cov_matrix), intent(in) :: C
!     REAL(cp), DIMENSION(:), INTENT(INOUT)  :: fb
!     REAL(cp), DIMENSION(:), INTENT(INOUT) :: gb
!     LOGICAL, OPTIONAL, INTENT(IN) :: die
!     !local variables
!     LOGICAL :: ll_die
!
!     IF( PRESENT(die) )THEN
!       ll_die = die
!     ELSE
!       ll_die = .TRUE.
!     END IF
!
!     SELECT CASE(C%im_nFac)
!       CASE (1)
!         CALL stop_program('apply_svd_square_covb is not yed defined for 1D')
!       CASE (2)
!         CALL cov_square_mvb(C%Ux, C%Dx, C%Uy, C%Dy, fb, gb)
!       CASE (default_nFact)!default value
!         IF(ll_die)THEN
!           CALL stop_program('In apply_svd_square_covb, covariance matrix is not loaded')
!         ELSE
!           fb=fb+gb
!           gb=0.0_cp
!         END IF
!     END SELECT
!   END SUBROUTINE apply_svd_square_covb
!
!   !> \brief apply the inverse of the square root of the covariance matrix to the vector f and result g
!   !! \param[in] C covariance matrix derived type
!   !! \param[in] f input vector
!   !! \param[in] g output vector
!   !! \param[in] [die] specify the behaviour of the program if the covariance matrix is not loaded. if TRUE, the program stops, otherwise g =f. The user should make sure that f and g have the same size.
!   !<
!   SUBROUTINE apply_svd_square_cov_inv(C, f, g, die)
!     type(svd_cov_matrix), intent(in) :: C
!     REAL(cp), DIMENSION(:), INTENT(IN)  :: f
!     REAL(cp), DIMENSION(:), INTENT(OUT) :: g
!     LOGICAL, OPTIONAL, INTENT(IN) :: die
!     !local variables
!     LOGICAL :: ll_die
!
!     IF( PRESENT(die) )THEN
!       ll_die = die
!     ELSE
!       ll_die = .TRUE.
!     END IF
!
!     SELECT CASE(C%im_nFac)
!       CASE (1)
!         CALL stop_program( 'apply_square_cov_inv is not yed defined for 1D' )
!       CASE (2)
!         CALL cov_square_inv_mv( C%Ux, C%Dx, C%Uy, C%Dy, f, g )
!       CASE ( default_nFact )!default value
!         IF (ll_die) THEN
!           CALL stop_program('In apply_svd_square_cov_inv, covariance matrix is not loaded')
!         ELSE
!           g=f
!         END IF
!     END SELECT
!   END SUBROUTINE apply_svd_square_cov_inv
!
! !end of public procedures
!
! !privates procedures
!   !> \brief Checks the parameters for the application of the covariance matrix or its inverse
!   !! @see cov_mv_2Fac for the description of the parameters
!   !<
!   SUBROUTINE cov_mv_check_2Fac(Ux, Dx, Uy, Dy, f, g)
!     REAL(cp), DIMENSION(:,:), INTENT(IN) :: Ux, Uy
!     REAL(cp), DIMENSION(:), INTENT(IN)  :: Dx, Dy, f
!     REAL(cp), DIMENSION(:), INTENT(IN) :: g
!     !local variables
!     INTEGER :: nx, ny, kx, ky, nf_predict
!
!     nx = SIZE(Ux,1)
!     ny = SIZE(Uy,1)
!     kx = SIZE(Ux,2)
!     ky = SIZE(Uy,2)
!     nf_predict = nx*ny
!     IF( (SIZE(Dx)/=kx).OR.(SIZE(Dy)/=ky).OR.(SIZE(f)/=nf_predict).OR.(SIZE(g)/=nf_predict) )THEN
!       CALL debug('Expecting data with following specification, based on the shape of Ux and Uy', tag=dALLWAYS)
!       CALL debug(kx, 'Dx of size :', tag=dALLWAYS)
!       CALL debug(ky, 'Dy of size :', tag=dALLWAYS)
!       CALL debug(nf_predict, 'f and g of size :', tag=dALLWAYS)
!       CALL debug('Get the following, check and try aigain', tag=dALLWAYS)
!       CALL debug(SIZE(Dx), 'SIZE(Dx)', tag=dALLWAYS)
!       CALL debug(SIZE(Dy), 'SIZE(Dy)', tag=dALLWAYS)
!       CALL debug(SIZE(f), 'SIZE(f)', tag=dALLWAYS)
!       CALL debug(SIZE(g), 'SIZE(g)', tag=dALLWAYS)
!       CALL debug('For debuging purpose, here are the shapes of Ux and Uy', tag=dALLWAYS)
!       CALL debug(SHAPE(Ux), 'SHAPE(Ux)', tag=dALLWAYS)
!       CALL debug(SHAPE(Uy), 'SHAPE(Uy)', tag=dALLWAYS)
!       CALL stop_program('Non consistant input for the application of the covariance matrix or its inverse')
!     END IF
!   END SUBROUTINE cov_mv_check_2Fac
!
!   !> \brief multiply a vector by the covariance matrix given as a Kronecker product of two covariance Bx and By. Each one being given by its SVD. Best suited for two physical dimensions problems where the covariance factorizes.
!   !! param[in] Ux matrix whose columns are the leading eigenvalues vector of the covariance matrix Bx in the direction x, size nx \times kx
!   !! param[in] Dx vector of leading eigenvalues in the direction x, size kx
!   !! param[in] Uy matrix whose columns are the leading eigenvalues vector of the covariance matrix By in the direction y, size ny \times ky
!   !! param[in] Dy vector of leading eigenvalues in the direction y, size ky
!   !! param[in] f input vector of size nx*ny (size of the discretization grid)
!   !! param[out] g output vector of size nx*ny
!   !! \details
!   !! computes g = Bf
!   !! this routine assumes either :
!   !! - the x-leading and and row major vectorization or
!   !! - the y-leading and the column major vectorization
!   !<
!   SUBROUTINE cov_mv_2Fac(Ux, Dx, Uy, Dy, f, g)
!     REAL(cp), DIMENSION(:,:), INTENT(IN) :: Ux, Uy
!     REAL(cp), DIMENSION(:), INTENT(IN)  :: Dx, Dy, f
!     REAL(cp), DIMENSION(:), INTENT(OUT) :: g
!     !local variables
!     REAL(cp), DIMENSION( SIZE(Dx)*SIZE(Dy) ) :: t1, t2
!
!     !checking the input parameters
!     CALL cov_mv_check(Ux, Dx, Uy, Dy, f, g)
!
!     CALL kron_matvec(Ux, Uy, f , t1, trans=.TRUE.)
!     CALL kron_matvec(Dx, Dy, t1, t2)
!     CALL kron_matvec(Ux, Uy, t2, g , trans=.FALSE.)
!   END SUBROUTINE cov_mv_2Fac
!
!   !> \brief Checks the parameters for the application of the covariance matrix or its inverse
!   !! @see cov_mv_2Fac for the description of the parameters
!   !<
!   SUBROUTINE cov_mv_check_default(U, D, f, g)
!     REAL(cp), DIMENSION(:,:), INTENT(IN) :: U
!     REAL(cp), DIMENSION(:), INTENT(IN)  :: D, f
!     REAL(cp), DIMENSION(:), INTENT(OUT) :: g
!     !local variables
!     INTEGER :: n, k
!
!     n = SIZE(U,1)
!     k = SIZE(U,2)
!     IF( (SIZE(D)/=k).OR.(SIZE(f)/=n).OR.(SIZE(g)/=n) )THEN
!       CALL debug('cov_mv_check_default: Expecting data with following specification, based on the shape of U', tag=dALLWAYS)
!       CALL debug(k, 'D of size :', tag=dALLWAYS)
!       CALL debug(n, 'f and g of size :', tag=dALLWAYS)
!       CALL debug('Get the following, check and try aigain', tag=dALLWAYS)
!       CALL debug(SIZE(D), 'SIZE(D)', tag=dALLWAYS)
!       CALL debug(SIZE(f), 'SIZE(f)', tag=dALLWAYS)
!       CALL debug(SIZE(g), 'SIZE(g)', tag=dALLWAYS)
!       CALL debug('For debuging purpose, here are the shapes of Ux and Uy', tag=dALLWAYS)
!       CALL debug(SHAPE(U), 'SHAPE(U)', tag=dALLWAYS)
!       CALL stop_program('Non consistant input for the application of the covariance matrix or its inverse')
!     END IF
!   END SUBROUTINE cov_mv_check_default
!
!   !> \brief multiply a vector by the covariance matrix given as its SVD
!   !! param[in] U matrix whose columns are the leading eigenvectors of the covariance matrix B, size n \times k
!   !! param[in] D vector of leading eigenvalues, size k
!   !! param[in] f input vector of size n (size of the discretization grid)
!   !! param[out] g output vector of size n
!   !! \details
!   !! computes g = Bf
!   !<
!   SUBROUTINE cov_mv_default(U, D, f, g)
!     REAL(cp), DIMENSION(:,:), INTENT(IN) :: U
!     REAL(cp), DIMENSION(:), INTENT(IN)  :: D, f
!     REAL(cp), DIMENSION(:), INTENT(OUT) :: g
!     !local variables
!     REAL(cp), DIMENSION( SIZE(D) ) :: t1, t2
!
!     !checking the input parameters
!     CALL cov_mv_check( U, D, f, g )
!
!     t1 = MATMUL(f, U) !transpose
!     t2 = D*t1
!     g  = MATMUL(U, t2)
!   END SUBROUTINE cov_mv_default
!
!   !> \brief multiplies a vector by the inverse of the covariance matrix, two physical dimensions problem
!   !! param[in] Ux matrix whose columns are the leading eigenvalues vector of the covariance matrix Bx in the direction x, size nx \times kx
!   !! param[in] Dx vector of leading eigenvalues in the direction x, size kx
!   !! param[in] Uy matrix whose columns are the leading eigenvalues vector of the covariance matrix By in the direction y, size ny \times ky
!   !! param[in] Dy vector of leading eigenvalues in the direction y, size ky
!   !! param[in] f input vector of size nx*ny (size of the discretization grid)
!   !! param[out] g output vector of size nx*ny
!   !! \details
!   !! computes g = B^-1 f
!   !! this routine assumes either :
!   !! - the x-leading and and row major vectorization or
!   !! - the y-leading and the column major vectorization
!   !<
!   SUBROUTINE cov_inv_mv_2Fac(Ux, Dx, Uy, Dy, f, g)
!     REAL(cp), DIMENSION(:,:), INTENT(IN) :: Ux, Uy
!     REAL(cp), DIMENSION(:), INTENT(IN)  :: Dx, Dy, f
!     REAL(cp), DIMENSION(:), INTENT(OUT) :: g
!     !local variables
!     REAL(cp), DIMENSION(SIZE(Dx)) :: DxInv
!     REAL(cp), DIMENSION(SIZE(Dy)) :: DyInv
!     REAL(cp), DIMENSION( SIZE(Dx)*SIZE(Dy) ) :: t1, t2
!
!     !checking the input parameters
!     CALL cov_mv_check(Ux, Dx, Uy, Dy, f, g)
!
!     DxInv = 1.0/Dx
!     DyInv = 1.0/Dy
!     CALL kron_matvec(Ux, Uy, f , t1, trans=.TRUE.)
!     CALL kron_matvec(DxInv, DyInv, t1, t2)
!     CALL kron_matvec(Ux, Uy, t2, g , trans=.FALSE.)
!   END SUBROUTINE cov_inv_mv_2Fac
!
!   !> \brief multiplies a vector by the inverse of the covariance matrix given as its SVD
!   !! param[in] U matrix whose columns are the leading eigenvalues vector of the covariance matrix B, size n \times k
!   !! param[in] D vector of leading eigenvalues, size k
!   !! param[in] f input vector of size n (size of the discretization grid)
!   !! param[out] g output vector of size n
!   !! \details
!   !! computes g = B^-1 f
!   !! this routine assumes either :
!   !! - the x-leading and and row major vectorization or
!   !! - the y-leading and the column major vectorization
!   !<
!   SUBROUTINE cov_inv_mv_default(U, D, f, g)
!     REAL(cp), DIMENSION(:,:), INTENT(IN) :: U
!     REAL(cp), DIMENSION(:), INTENT(IN)  :: D, f
!     REAL(cp), DIMENSION(:), INTENT(OUT) :: g
!     !local variables
!     REAL(cp), DIMENSION(SIZE(D)) :: t1, t2
!
!     !checking the input parameters
!     CALL cov_mv_check( U, D, f, g )
!
!     t1 = MATMUL(f, U) !transpose
!     t2 = t1/D
!     g  = MATMUL(U, t2)
!   END SUBROUTINE cov_inv_mv_default
!
!   !> \brief Checks the parameters for the application of the square root of the covariance matrix or its inverse
!   !! @see cov_square_mv_default for the description of the parameters
!   !<
!   SUBROUTINE cov_square_mv_check_default(U, D, f, g, inv)
!     REAL(cp), DIMENSION(:,:), INTENT(IN) :: U
!     REAL(cp), DIMENSION(:), INTENT(IN)  :: D, f
!     REAL(cp), DIMENSION(:), INTENT(OUT) :: g
!     LOGICAL, OPTIONAL :: inv
!     !local variables
!     INTEGER :: n, k, nf_predict, ng_predict
!     LOGICAL :: ll_inv
!
!     IF( PRESENT(inv) )THEN
!       ll_inv = inv
!     ELSE
!       ll_inv = .FALSE.
!     END IF
!     n = SIZE(U,1)
!     k = SIZE(U,2)
!     IF(ll_inv)THEN
!       nf_predict = n
!       ng_predict = k
!     ELSE
!       nf_predict = k
!       ng_predict = n
!     END IF
!     IF( (SIZE(D)/=k).OR.(SIZE(f)/=nf_predict).OR.(SIZE(g)/=ng_predict) )THEN
!       CALL debug('Expecting data with following specification, based on the shape of Ux and Uy', tag=dALLWAYS)
!       CALL debug(k, 'D of size :', tag=dALLWAYS)
!       CALL debug(nf_predict, 'f of size :', tag=dALLWAYS)
!       CALL debug(ng_predict, 'g of size :', tag=dALLWAYS)
!       CALL debug('Get the following, check and try aigain', tag=dALLWAYS)
!       CALL debug(SIZE(D), 'SIZE(Dx)', tag=dALLWAYS)
!       CALL debug(SIZE(f), 'SIZE(f)', tag=dALLWAYS)
!       CALL debug(SIZE(g), 'SIZE(g)', tag=dALLWAYS)
!       CALL debug('For debuging purpose, here are the shapes of Ux and Uy', tag=dALLWAYS)
!       CALL debug(SHAPE(U), 'SHAPE(U)', tag=dALLWAYS)
!       IF(ll_inv) CALL debug('The purpose is to compute apply the inverse square root', tag=dALLWAYS)
!       CALL stop_program('Non consistant input for the application of the covariance matrix or its inverse')
!     END IF
!   END SUBROUTINE cov_square_mv_check_default
!
!   !> \brief Checks the parameters for the application of the square root of the covariance matrix or its inverse
!   !! @see cov_square_mv_2Fac for the description of the parameters
!   !<
!   SUBROUTINE cov_square_mv_check_2Fac(Ux, Dx, Uy, Dy, f, g, inv)
!     REAL(cp), DIMENSION(:,:), INTENT(IN) :: Ux, Uy
!     REAL(cp), DIMENSION(:), INTENT(IN)  :: Dx, Dy, f
!     REAL(cp), DIMENSION(:), INTENT(IN) :: g
!     LOGICAL, OPTIONAL :: inv
!     !local variables
!     INTEGER :: nx, ny, kx, ky, nf_predict, ng_predict
!     LOGICAL :: ll_inv
!
!     IF( PRESENT(inv) )THEN
!       ll_inv = inv
!     ELSE
!       ll_inv = .FALSE.
!     END IF
!     nx = SIZE(Ux,1)
!     ny = SIZE(Uy,1)
!     kx = SIZE(Ux,2)
!     ky = SIZE(Uy,2)
!     IF(ll_inv)THEN
!       nf_predict = nx*ny
!       ng_predict = kx*ky
!     ELSE
!       nf_predict = kx*ky
!       ng_predict = nx*ny
!     END IF
!     IF( (SIZE(Dx)/=kx).OR.(SIZE(Dy)/=ky).OR.(SIZE(f)/=nf_predict).OR.(SIZE(g)/=ng_predict) )THEN
!       CALL debug('Expecting data with following specification, based on the shape of Ux and Uy', tag=dALLWAYS)
!       CALL debug(kx, 'Dx of size :', tag=dALLWAYS)
!       CALL debug(ky, 'Dy of size :', tag=dALLWAYS)
!       CALL debug(nf_predict, 'f of size :', tag=dALLWAYS)
!       CALL debug(ng_predict, 'g of size :', tag=dALLWAYS)
!       CALL debug('Get the following, check and try aigain', tag=dALLWAYS)
!       CALL debug(SIZE(Dx), 'SIZE(Dx)', tag=dALLWAYS)
!       CALL debug(SIZE(Dy), 'SIZE(Dy)', tag=dALLWAYS)
!       CALL debug(SIZE(f), 'SIZE(f)', tag=dALLWAYS)
!       CALL debug(SIZE(g), 'SIZE(g)', tag=dALLWAYS)
!       CALL debug('For debuging purpose, here are the shapes of Ux and Uy', tag=dALLWAYS)
!       CALL debug(SHAPE(Ux), 'SHAPE(Ux)', tag=dALLWAYS)
!       CALL debug(SHAPE(Uy), 'SHAPE(Uy)', tag=dALLWAYS)
!       IF(ll_inv) CALL debug('The purpose is to compute apply the inverse square root', tag=dALLWAYS)
!       CALL stop_program('Non consistant input for the application of the covariance matrix or its inverse')
!     END IF
!   END SUBROUTINE cov_square_mv_check_2Fac
!
!   !> \brief multiplies a vector by the square root of the covariance matrix given as its SVD
!   !! param[in] U matrix whose columns are the leading eigenvalues vector of the covariance matrix B, size n \times k
!   !! param[in] D vector of leading eigenvalues, size k
!   !! param[in] f input vector of size n (size of the discretization grid)
!   !! param[out] g output vector of size n
!   !! \details
!   !! computes g = B^0.5 f
!   !<
!   SUBROUTINE cov_square_mv_default(U, D, f, g)
!     REAL(cp), DIMENSION(:,:), INTENT(IN) :: U
!     REAL(cp), DIMENSION(:), INTENT(IN)  :: D, f
!     REAL(cp), DIMENSION(:), INTENT(OUT) :: g
!     !local variables
!     REAL(cp), DIMENSION( SIZE(D) ) :: t
!
!     !checking the input parameters
!     CALL cov_square_mv_check(U, D, f, g, inv=.FALSE.)
!
!     t = SQRT(D)*f
!     g = MATMUL(U, t)
!   END SUBROUTINE cov_square_mv_default
!
!   !> \brief multiplies a vector by the square root of the covariance matrix, two physical dimensions problem
!   !! param[in] Ux matrix whose columns are the leading eigenvalues vector of the covariance matrix Bx in the direction x, size nx \times kx
!   !! param[in] Dx vector of leading eigenvalues in the direction x, size kx
!   !! param[in] Uy matrix whose columns are the leading eigenvalues vector of the covariance matrix By in the direction y, size ny \times ky
!   !! param[in] Dy vector of leading eigenvalues in the direction y, size ky
!   !! param[in] f input vector of size nx*ny (size of the discretization grid)
!   !! param[out] g output vector of size nx*ny
!   !! \details
!   !! computes g = B^0.5 f
!   !! this routine assumes either :
!   !! - the x-leading and and row major vectorization or
!   !! - the y-leading and the column major vectorization
!   !<
!   SUBROUTINE cov_square_mv_2Fac(Ux, Dx, Uy, Dy, f, g)
!     REAL(cp), DIMENSION(:,:), INTENT(IN) :: Ux, Uy
!     REAL(cp), DIMENSION(:), INTENT(IN)  :: Dx, Dy, f
!     REAL(cp), DIMENSION(:), INTENT(OUT) :: g
!     !local variables
!     REAL(cp), DIMENSION(SIZE(Dx)) :: Dx05
!     REAL(cp), DIMENSION(SIZE(Dy)) :: Dy05
!     REAL(cp), DIMENSION( SIZE(Dx)*SIZE(Dy) ) :: t
!
!     !checking the input parameters
!     CALL cov_square_mv_check_2Fac(Ux, Dx, Uy, Dy, f, g, inv=.FALSE.)
!
!     Dx05 = SQRT(Dx)
!     Dy05 = SQRT(Dy)
!
!     CALL kron_matvec(Dx05, Dy05, f, t)
!     CALL kron_matvec(Ux, Uy, t, g , trans=.FALSE.)
!   END SUBROUTINE cov_square_mv_2Fac
!
!   !> adjoint subroutine, used mathmatical formulation of the adjoint with linear operators
!   SUBROUTINE cov_square_mvb_default(U, D, fb, gb)
!     REAL(cp), DIMENSION(:,:), INTENT(IN) :: U
!     REAL(cp), DIMENSION(:), INTENT(IN)  :: D
!     REAL(cp), DIMENSION(:), INTENT(INOUT) :: gb, fb
!     !local variables
!     REAL(cp), DIMENSION( SIZE(D) ) :: tb
!
!     !checking the input parameters
!     CALL cov_square_mv_check(U, D, fb, gb, inv=.FALSE.)
!
!     tb = tb + MATMUL(gb, U)
!     gb = 0.0_cp
!     fb = fb + SQRT(D)*tb
!     tb = 0.0_cp
!   END SUBROUTINE cov_square_mvb_default
!
!   !> adjoint subroutine, used mathmatical formulation of the adjoint with linear operators
!   SUBROUTINE cov_square_mvb_2Fac(Ux, Dx, Uy, Dy, fb, gb)
!     REAL(cp), DIMENSION(:,:), INTENT(IN) :: Ux, Uy
!     REAL(cp), DIMENSION(:), INTENT(IN)  :: Dx, Dy
!     REAL(cp), DIMENSION(:), INTENT(INOUT) :: gb, fb
!     !local variables
!     REAL(cp), DIMENSION(SIZE(Dx)) :: Dx05
!     REAL(cp), DIMENSION(SIZE(Dy)) :: Dy05
!     REAL(cp), DIMENSION( SIZE(Dx)*SIZE(Dy) ) :: tb
!     REAL(cp), DIMENSION( SIZE(fb) ) :: tmp
!
!     !checking the input parameters
!     CALL cov_square_mv_check_2Fac(Ux, Dx, Uy, Dy, fb, gb, inv=.FALSE.)
!     !initialization
!
!     Dx05 = SQRT(Dx)
!     Dy05 = SQRT(Dy)
!
!     CALL kron_matvec( Ux, Uy, gb, tb , trans=.TRUE. )!/!\optimised, knowing that tb is zero
!     gb = 0.0_cp
!     CALL kron_matvec( Dx05, Dy05, tb, tmp )! No need for transpose with diagonal square matrices
!     fb = fb + tmp
!     tb = 0.0_cp
!   END SUBROUTINE cov_square_mvb_2Fac
!
!   !> \brief multiplies a vector by the inverse of the square root of the covariance matrix, two physical dimensions problem
!   !! param[in] Ux matrix whose columns are the leading eigenvalues vector of the covariance matrix Bx in the direction x, size nx \times kx
!   !! param[in] Dx vector of leading eigenvalues in the direction x, size kx
!   !! param[in] Uy matrix whose columns are the leading eigenvalues vector of the covariance matrix By in the direction y, size ny \times ky
!   !! param[in] Dy vector of leading eigenvalues in the direction y, size ky
!   !! param[in] f input vector of size nx*ny (size of the discretization grid)
!   !! param[out] g output vector of size nx*ny
!   !! \details
!   !! computes g = B^-0.5 f
!   !! this routine assumes either :
!   !! - the x-leading and and row major vectorization or
!   !! - the y-leading and the column major vectorization
!   !<
!   SUBROUTINE cov_square_inv_mv_2Fac(Ux, Dx, Uy, Dy, f, g)
!     REAL(cp), DIMENSION(:,:), INTENT(IN) :: Ux, Uy
!     REAL(cp), DIMENSION(:), INTENT(IN)  :: Dx, Dy, f
!     REAL(cp), DIMENSION(:), INTENT(OUT) :: g
!     !local variables
!     REAL(cp), DIMENSION(SIZE(Dx)) :: DxInv
!     REAL(cp), DIMENSION(SIZE(Dy)) :: DyInv
!     REAL(cp), DIMENSION( SIZE(Dx)*SIZE(Dy) ) :: t
!
!     !checking the input parameters
!     CALL cov_square_mv_check_2Fac(Ux, Dx, Uy, Dy, f, g, inv=.TRUE.)
!
!     DxInv = SQRT(1.0/Dx)
!     DyInv = SQRT(1.0/Dy)
!     CALL kron_matvec(Ux, Uy, f, t , trans=.TRUE.)
!     CALL kron_matvec(DxInv, DyInv, t, g)
!   END SUBROUTINE cov_square_inv_mv_2Fac
!
!   !> \brief multiplies a vector by the inverse of the square root of the covariance matrix given as its SVD
!   !! param[in] U matrix whose columns are the leading eigenvalues vector of the covariance matrix B, size n \times k
!   !! param[in] D vector of leading eigenvalues, size k
!   !! param[in] f input vector of size n (size of the discretization grid)
!   !! param[out] g output vector of size n
!   !! \details
!   !! computes g = B^-0.5 f
!   !<
!   SUBROUTINE cov_square_inv_mv_default(U, D, f, g)
!     REAL(cp), DIMENSION(:,:), INTENT(IN) :: U
!     REAL(cp), DIMENSION(:), INTENT(IN)  :: D, f
!     REAL(cp), DIMENSION(:), INTENT(OUT) :: g
!     !local variables
!     REAL(cp), DIMENSION( SIZE(D) ) :: t
!
!     !checking the input parameters
!     CALL cov_square_mv_check(U, D, f, g, inv=.TRUE.)
!
!     t = MATMUL(f, U)
!     g = t/SQRT(D)
!   END SUBROUTINE cov_square_inv_mv_default
!
!   ! ***************************************************
!   ! **************** debugging routines ***************
!   ! ***************************************************
!
!   !> \brief Check the conformity of the shape of matrices forming an SVD
!   !! the program stops in case of non conformity
!   !! \param[in] nU size of each singular vector
!   !! \param[in] kU number singular vectors
!   !! \param[in] kD1 number of singular values or on size of the diagonal matrix
!   !! \param[in] kD2 optional (default 1) other size of the diagonal matrix
!   !! kD1 x kD2 is the shape of the diagonal matrix of singular values. It can be stored as a full matrix or as a vector.
!   !<
!   SUBROUTINE SVD_check(nU, kU, kD1, kD2)
!     INTEGER, INTENT(IN) :: nU, kU, kD1
!     INTEGER, OPTIONAL, INTENT(IN) :: kD2
!     !local variables
!     INTEGER :: il_kD2
!
!     IF( PRESENT(kD2)  )THEN
!       il_kD2 = kD2
!     ELSE
!       il_kD2 = 1
!     END IF
!
!     IF( ( (kD1 == il_kD2).OR.( MIN(kD1, il_kD2) == 1 ) ) &
!       .AND.( MAX(kD1, il_kD2) == kU ).AND.( nU >= kU ) )THEN
!       !good configuration, nothing to do
!     ELSE
!       !bad configuration, message and stop the program
!       CALL debug( 'Expecting SVD UDU^t to have the following specifications', tag=dALLWAYS )
!       CALL debug( ' - U: number of rows greater or equal to the number of columns', tag=dALLWAYS )
!       CALL debug( ' - D: vector or square matrix', tag=dALLWAYS )
!       CALL debug( ' - MIN size of U equal to the MAX size of D', tag=dALLWAYS )
!       CALL debug( 'For debuging purpose, here are the shapes of U and D', tag=dALLWAYS )
!       CALL debug( (/nU, kU/), ' - SHAPE(U) = ', tag=dALLWAYS )
!       CALL debug( (/kD1, il_kD2/), ' - SHAPE(D) = ', tag=dALLWAYS )
!       CALL stop_program( 'Non consistant SVD' )
!     END IF
!   END SUBROUTINE SVD_check
! !end of privates procedures
! END MODULE covariance

max_stars_repo_name: souopgui/da-toolbox

max_issues_repo_name: souopgui/da-toolbox

max_forks_repo_name: souopgui/da-toolbox


text:
	SUBROUTINE  DM_GETF ( iflno, iret )
C************************************************************************
C* DM_GETF								*
C*									*
C* This subroutine gets the next file number to use.			*
C*									*
C* DM_GETF  ( IFLNO, IRET )						*
C*									*
C* Output parameters:							*
C*	IFLNO		INTEGER		File number			*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -3 = too many files open	*
C**									*
C* Log:									*
C* M. desJardins/GSFC	6/86						*
C* T. Piper/GSC		11/98	Updated prolog				*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
	INCLUDE		'GMBDTA.CMN'
	INCLUDE		'dmcmn.cmn'
C-----------------------------------------------------------------------
C*	Initialize and find first available file number.
C
	iflno = 0
	iret  = 0
	DO  i = MMFILE, 1, -1
	    IF  ( lundm (i) .lt. 0 ) iflno = i
	END DO
	IF  ( iflno .eq. 0 ) iret = -3
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: sgdecker/pygempak


text:
module mod_fillps
  use mod_types
  implicit none
  private
  public fillps
  contains
  subroutine fillps(n,dli,dzfi,dti,up,vp,wp,p)
    !@cuf use cudafor
    !
    !  fill the right-hand side of the Poisson equation for the correction pressure.
    !
    !  the discrete divergence is:
    !
    !  w(i,j,k)-w(i,j,k-1)   v(i,j,k)-v(i,j-1,k)   u(i,j,k)-u(i-1,j,k)
    !  ------------------- + ------------------- + -------------------  = div
    !          dz                    dy                    dx
    !
    implicit none
    integer , intent(in), dimension(3) :: n
    real(rp), intent(in), dimension(3) :: dli
    real(rp), intent(in), dimension(0:) :: dzfi
    real(rp), intent(in) :: dti
    real(rp), intent(in ), dimension(0:,0:,0:) :: up,vp,wp
    real(rp), intent(out), dimension(0:,0:,0:) :: p
    real(rp) :: dtidxi,dtidyi!,dtidzi
    !real(rp), dimension(0:n(3)+1) :: dtidzfi
#ifdef USE_CUDA
    attributes(managed) :: up,vp,wp,p,dzfi!,dtidzfi
    integer :: istat
#endif
    integer :: i,j,k,im,jm,km
    !
    dtidxi = dti*dli(1)
    dtidyi = dti*dli(2)
#ifdef USE_CUDA
   !$cuf kernel do(3) <<<*,*>>>
#else
    !$OMP PARALLEL DO DEFAULT(none) &
    !$OMP SHARED(n,p,up,vp,wp,dti,dzfi,dtidyi,dtidxi) &
    !$OMP PRIVATE(i,j,k,im,jm,km)
#endif
    do k=1,n(3)
      km = k-1
      do j=1,n(2)
        jm = j-1
        do i=1,n(1)
          im = i-1
          p(i,j,k) = ( &
                      (wp(i,j,k)-wp(i,j,km))*dti*dzfi(k)+ &
                      (vp(i,j,k)-vp(i,jm,k))*dtidyi    + &
                      (up(i,j,k)-up(im,j,k))*dtidxi    )
        enddo
      enddo
    enddo
#ifndef USE_CUDA
    !$OMP END PARALLEL DO
#endif
    !@cuf istat=cudaDeviceSynchronize()
    !
    return
  end subroutine fillps
end module mod_fillps

max_stars_repo_name: maxcuda/CaNS

max_issues_repo_name: M-Abbaszadeh/CaNS

max_forks_repo_name: M-Abbaszadeh/CaNS


text:
!*robodoc*u* tests/testABK
!  NAME
!    testABK
!  SYNOPSIS

!$Id: testABK.f90 380 2017-03-22 11:03:09Z mexas $

program testABK

!  PURPOSE
!    Checking: cgca_av, cgca_dv, cgca_gv
!  DESCRIPTION
!    Checking and timing grain volume array calculation on Cray.
!    Use with CRAY only!
!  NOTES
!    The program must be called with 2 command line arguments,
!    both positive integers. These are codimensions along 1 and 2.
!    The number of images must be such that
!    codimension3 = num_images()/( codimension1 * codimension3 )
!    is a positive integer. Example:
!      ./testABK.x 2 2
!    which will make the third codimension equal to 16/(2*2)=4.
!  AUTHOR
!    Anton Shterenlikht
!  COPYRIGHT
!    See LICENSE
!  USES
!    cgca testaux
!  USED BY
!    Part of CGPACK test suite
!  SOURCE

use testaux

implicit none

 logical(kind=ldef),parameter :: yesdebug=.true., nodebug=.false., &
   periodicbc=.true.
 real,parameter :: gigabyte=real(2**30), resolution=1.0e-5

 integer(kind=idef) :: l1,u1,l2,u2,l3,u3,col1,cou1,col2,cou2,col3,cou3, &
   nuc,    & ! number of nuclei in the model
   nimages, img, codim(3)[*], i
 integer(kind=iarr),allocatable :: space(:,:,:,:)[:,:,:]
 integer(kind=ilrg),allocatable :: grainvol(:)[:,:,:]
 integer(kind=ilrg) :: icells, mcells

 logical(kind=ldef) :: solid

 real :: image_storage

 integer :: pat_status

!*********************************************************************72
! first executable statement

nimages = num_images()
img = this_image()

! do a check on image 1
if ( img .eq. 1 ) then
 call getcodim(nimages,codim)
 ! print a banner
 call banner("ABK")
 ! print the parameter values
 call cgca_pdmp
 write (*,'(a,i0,a)') "running on ", nimages, " images in a 3D grid"
end if

! all images read codim from image 1
if ( img .eq. 1 ) then
 sync images(*)
else
 sync images(1)
 codim(:) = codim(:)[1]
end if

l1=1
l2=l1
l3=l1

! The array size is only controlled by this value
! in this program.
u1=128
u2=u1
u3=u1

col1=1
cou1=codim(1)-col1+1
col2=1
cou2=codim(2)-col2+1
col3=1
cou3=codim(3)-col3+1

! total number of cells in a coarray
icells = int(u1-l1+1,kind=ilrg) * int(u2-l2+1,kind=ilrg) * &
  int(u3-l3+1,kind=ilrg)

! total number of cells in the model
mcells = icells * int(codim(1),kind=ilrg) * int(codim(2),kind=ilrg) * &
  int(codim(3),kind=ilrg)

! total number of nuclei
nuc = int( resolution * real(mcells) )

if ( img .eq. 1 ) then
  write (*,'(a,2(i0,":",i0,","),i0,":",i0,")")') &
    "bounds: (",l1,u1,l2,u2,l3,u3
  write (*,'(a,2(i0,":",i0,","),i0,":",i0,")")') &
    "cobounds: (",col1,cou1,col2,cou2,col3,cou3

  ! An absolute minimum of storage, in GB, per image.
  ! A factor of 2 is used because will call _sld, which
  ! allocates another array of the same size and kind as
  ! coarray.
  image_storage = real(2 * icells*storage_size(space)/8)/gigabyte

  write (*,'(a,i0,a)') "Each image has ",icells, " cells"
  write (*,'(a,i0,a)') "The model has ", mcells, " cells"
  write (*,'(a,i0,a)') "The model has ", nuc, " nuclei"
  write (*,'(a,es9.2,a)') "Each image will use at least ", &
    image_storage, " GB memory"
end if

! initialise random number seed
call cgca_irs(nodebug)

! allocate space and grain volume coarrays
call cgca_as(l1,u1,l2,u2,l3,u3,col1,cou1,col2,cou2,col3,1,space)
call cgca_av(1,nuc,col1,cou1,col2,cou2,col3,grainvol)

! initialise coarray
space = cgca_liquid_state

! populate nuclei
call cgca_nr(space,nuc,nodebug)

! solidify
call cgca_sld(space,periodicbc,0,1,solid)

! calculate volumes with my routine
call pat_region_begin(1,"_gv",pat_status)
call cgca_gv(space,grainvol)
call pat_region_end(1,pat_status)

! dump grain volumes
if ( img .eq. 1 ) then
 write (*,*) "results from cgca_gv"
 do i=lbound(grainvol,dim=1),ubound(grainvol,dim=1)
  write (*,"(i0,tr1,i0)") i, grainvol(i)
 end do
end if

! now calculate volumes using CO_SUM intrinsic
call pat_region_begin(2,"_gvl+co_sum",pat_status)
call cgca_gvl(space,grainvol)
call co_sum(grainvol)
call pat_region_end(2,pat_status)

! dump grain volumes
if ( img .eq. 1 ) then
 write (*,*) "results from cgca_gvl + co_sum"
 do i=lbound(grainvol,dim=1),ubound(grainvol,dim=1)
  write (*,"(i0,tr1,i0)") i, grainvol(i)
 end do
end if

! deallocate all arrays
call cgca_ds(space)
call cgca_dv(grainvol)

end program testABK

!*roboend*

max_stars_repo_name: lcebaman/casup

max_issues_repo_name: lcebaman/casup

max_forks_repo_name: lcebaman/casup


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  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, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

!
! Character function to concat two strings into a new string.
! The arguments of this function are similar to the EXEC2 function
! &CONCAT OF.
!
      CHARACTER*(*) FUNCTION CONCAT (WORD1,WORD2,MAXLE1,MAXLE2)
      CHARACTER*(*) WORD1,WORD2
      character*1   DUMMY
      INTEGER MAXLE1, MAXLE2, LENSTR, STOSH
!
! Obtain the lengths of WORD1 and WORD2.
! Use STOSH to write all significant characters into CONCAT.
!
      CONCAT = DUMMY
      LEN1 = LENSTR (WORD1,MAXLE1)
      LEN2 = LENSTR (WORD2,MAXLE2)
      IRC = STOSH(WORD1, 1, LEN1, CONCAT, 1, LENOUT)
      IRC = STOSH(WORD2, 1, LEN2, CONCAT, LENOUT+1, LENTOT)
      RETURN
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
        SUBROUTINE BCKREC ( FILE )
      INCLUDE 'DSIOF.COM'
        INTEGER FILE
        NAME = FILE
        CALL DSGEFL
        CALL DSBRC1
        CALL DSSDCB
        RETURN
        END


text: x <= a or (a and b);


text:
!=======================================================================
!  SALUS SOIL EVAPORATION MODULE - File ESR_SoilEvap.for
!=======================================================================
!  ESR_SoilEvap, Subroutine, J. Ritchie, C. Porter
!
!  Calculates actual soil evaporation (ES, mm/d) based on method
!  described in:
!
!  Ritchie, J.T., C.H. Porter, J.Judge, J.W.Jones, A.A. Suleiman. 2009.
!    Application of a functional model for simulation of soil evaporation
!    and water redistribution.  Division S-1 -- soil Physics;
!    Soil Science Society of America. in review.
!
!  and
!
!  Suleiman, A.A., J.T.Ritchie. Modeling Soil Water Redistribution
!    during Second-Stage Evaporation. Division S-1 -- soil Physics;
!    Soil Science Society of America. Vol. 67, No. 2. March-apr 2003.
!
!  This routine takes the place of SOILEV and UPFLOW.
!-----------------------------------------------------------------------
!  REVISION HISTORY
!  05/03/2005 JTR/CHP Written
!  07/05/2006 JTR/CHP Modified calculation of C_dry, function of depth
!                         only.  Eliminate A and B coefficients.
!  04/08/2008 JTR modification for wet profile
!  05/29/2008 JTR added intermediate profile case
!  10/02/2008 CHP/JTR changed depth for determining evaporation case
!                     from 50 cm to 100 cm.
!-----------------------------------------------------------------------
!  Called by: SPAM
!=======================================================================
      SUBROUTINE ESR_SoilEvap(
     &   EOS, SOILPROP, SW, SWDELTS,                      !Input
     &   ES, ES_LYR, SWDELTU, UPFLOW)                     !Output

!-----------------------------------------------------------------------
      USE ModuleDefs; USE ModuleData
      IMPLICIT NONE
      SAVE

!     ------------------------------------------------
!     Interface Variables:
      REAL, INTENT(IN) :: EOS          !Potential soil evap (mm/d)
      REAL, INTENT(IN) :: SW(NL)       !Soil water content (cm3/cm3)
      REAL, INTENT(IN) :: SWDELTS(NL)  !Rate of drainage (cm3/cm3)
      TYPE (SoilType), INTENT(IN) :: SOILPROP !Soil properties

      REAL, INTENT(OUT):: ES           !Actual soil evaporation (mm/d)
      REAL, INTENT(OUT):: SWDELTU(NL)  !Change in soil water (cm3/cm3)
      REAL, INTENT(OUT):: UPFLOW(NL)   !Flow or N transport (cm/d)
      REAL, INTENT(OUT):: ES_LYR(NL)   !Actual soil evap by layer (mm/d)
!     UPFLOW(1:NL) refers to water which moves up from layer L to
!       layer L-1, and includes upflow from lower layers.
!     ------------------------------------------------

!      CHARACTER*12, PARAMETER :: ERRKEY = 'SAL_SoilEvap'
      INTEGER L, NLAYR, ProfileType
      REAL A, B, RedFac, SW_threshold
      REAL, DIMENSION(NL) :: DLAYR, DS, DUL, LL, MEANDEP
      REAL, DIMENSION(NL) :: SWAD, SWTEMP, SW_AVAIL, ES_Coef

!-----------------------------------------------------------------------
!     ProfileType:
!     1 = Wet: SW > DUL in at least one layer in top 100 cm and
!           SW > SW_threshold in top layer
!     2 = Intermediate: wet, but SW < SW_threshold in top layer
!     3 = Dry: SW < DUL in all layers in top 100 cm
!-----------------------------------------------------------------------

      DLAYR = SOILPROP % DLAYR
      DS    = SOILPROP % DS
      DUL   = SOILPROP % DUL
      LL    = SOILPROP % LL
      NLAYR = SOILPROP % NLAYR

      ES = 0.0

!**********************************************************************
!     NEW 4/18/2008
      ProfileType = 3   !assume dry profile until proven wet
      DO L = 1, NLAYR
!       Air dry water content
        SWAD(L) = 0.30 * LL(L) !JTR 11/28/2006

!       Mean depth for each soil layer
        MEANDEP(L) = DS(L) - DLAYR(L) / 2.               !cm

!       Pseudo-integraton step
!       If increase in SW due to rain or irrigation, include half
        IF (SWDELTS(L) > 0.0) THEN
          SWTEMP(L) = SW(L) + 0.5 * SWDELTS(L)
        ELSE
!         If decrease in SW due to drainage, include all
          SWTEMP(L) = SW(L) + SWDELTS(L)
        ENDIF

!       If any layer in top 100 cm is wet, use wet profile method
        IF (MEANDEP(L) < 100. .AND. SWTEMP(L) > DUL(L)) THEN
          ProfileType = 1
        ENDIF
      ENDDO

!     If wet profile, check for top layer SW below threshold.
      IF (ProfileType == 1) THEN
!       SW_threshold = DUL(1) - 0.05 !/ 0.13 * (DUL(1) - LL(1))
!       JTR 6/4/2008
!       Threshold WC = 0.275*DUL +1.165*DUL^2 + (1.2*DUL^3.75)*depth (center)
        SW_threshold = 0.275*DUL(1) + 1.165*DUL(1)*DUL(1) +
     &          (1.2*DUL(1)**3.75)*MEANDEP(1)
!       chp 6/4/2008 use DUL - 0.05, like before, but limit to air dry
!        SW_threshold = MAX(SWAD(1), DUL(1) - 0.05)
        IF (SWTEMP(1) < SW_threshold) THEN
          ProfileType = 2
        ENDIF
      ENDIF

      DO L = 1, NLAYR
!-----------------------------------------------------------------------
        SELECT CASE (ProfileType)

!       Dry profile
        CASE (3)
!         Depth-dependant coefficients based on Ritchie spreadsheet 11/29/2006
          A =  0.5  + 0.24 * DUL(L)
          B = -2.04 + 0.20 * DUL(L)
          ES_Coef(L) = A * MEANDEP(L) ** B

!       Equilibrium profile
        CASE (2)
          ES_Coef(L) = 0.011   !for all depths

!       Wet profile
        CASE (1)
!         Ritchie spreadsheet of 5/28/08
          A = 0.26  !6/20/08  A = 0.14  !6/2/08  A = 0.42   !4/18/08
          B = -0.70 !6/20/08  B = -0.46 !6/2/08  B = -0.73  !4/18/08
          ES_Coef(L) = A * MEANDEP(L) ** B !function, no integration

        END SELECT
!-----------------------------------------------------------------------

        SWDELTU(L) = -(SWTEMP(L) - SWAD(L)) * ES_Coef(L) !mm3/mm3

!       Limit to available water
        SW_AVAIL(L) = SW(L) + SWDELTS(L) - SWAD(L)
        IF (-SWDELTU(L) > SW_AVAIL(L)) THEN
          SWDELTU(L) = -SW_AVAIL(L)                   !mm3/mm3
        ENDIF

!       Limit to negative values (decrease SW)
        SWDELTU(L) = AMIN1(0.0, SWDELTU(L))

!       Aggregate soil evaporation from each layer
        ES_LYR(L) = -SWDELTU(L) * DLAYR(L) * 10.      !mm
        ES = ES + ES_LYR(L)                           !profile sum (mm)
      ENDDO

!     Limit total profile soil evaporation to potential soil evaporation
      RedFac = 1.0
      If (ES > EOS) Then
        RedFac = EOS / ES
        ES_LYR = ES_LYR * RedFac
        SWDELTU = SWDELTU * RedFac
        ES = EOS
      End If

      UPFLOW = 0.0
      UPFLOW(NLAYR) = ES_LYR(NLAYR) / 10.
      DO L = NLAYR-1, 1, -1
        UPFLOW(L) = UPFLOW(L+1) + ES_LYR(L) / 10.     !cm/d
      ENDDO

!-----------------------------------------------------------------------
      RETURN
      END SUBROUTINE ESR_SoilEvap
!=======================================================================

ext: for

max_stars_repo_name: vstocca/dssat-csm-os

max_issues_repo_name: vstocca/dssat-csm-os

max_forks_repo_name: vstocca/dssat-csm-os


text:
module median_stats_mod
    use :: quicksort_mod
    implicit none

    private

    type, public :: descriptive_stats_t
        private
        integer :: nr_values = 0
        real :: sum = 0.0, sum2 = 0.0
        real, dimension(:), allocatable :: values
    end type descriptive_stats_t

    public :: init_stats, add_value, get_nr_values, get_mean, get_stddev, &
        get_median, finalize_stats

contains

    subroutine init_stats(stats, nr_values)
        use, intrinsic :: iso_fortran_env, only : error_unit
        type(descriptive_stats_t), intent(inout) :: stats
        integer, value :: nr_values
        integer :: status

        allocate(stats%values(nr_values), stat=status)
        if (status /= 0) then
            write (unit=error_unit, fmt='(A, I0, A)') &
                'error: can not allocated memory for ', nr_values, ' values'
            stop 3
        end if
    end subroutine init_stats

    subroutine add_value(stats, new_value)
        use, intrinsic :: iso_fortran_env, only : error_unit
        implicit none
        type(descriptive_stats_t), intent(inout) :: stats
        real, value :: new_value

        stats%nr_values = stats%nr_values + 1
        if (stats%nr_values > size(stats%values)) then
            write (unit=error_unit, fmt='(A, I0, A)') &
                'error: capacity exceeded, maximum ', size(stats%values), &
                ' can be processed'
            stop 4
        end if
        stats%sum = stats%sum + new_value
        stats%sum2 = stats%sum2 + new_value**2
        stats%values(stats%nr_values) = new_value
    end subroutine add_value

    function get_nr_values(stats) result(nr_values)
        implicit none
        type(descriptive_stats_t), intent(in) :: stats
        integer :: nr_values

        nr_values = stats%nr_values
    end function get_nr_values

    function get_mean(stats) result(mean_value)
        use, intrinsic :: iso_fortran_env, only : error_unit
        implicit none
        type(descriptive_stats_t), intent(inout) :: stats
        real :: mean_value

        if (stats%nr_values == 0) then
            write (unit=error_unit, fmt='(A)') &
                'error: too few data value to compute mean'
            stop 1
        end if
        associate (n => stats%nr_values, sum => stats%sum)
            mean_value = sum/n
        end associate
    end function get_mean

    function get_stddev(stats) result(stddev_value)
        use, intrinsic :: iso_fortran_env, only : error_unit
        implicit none
        type(descriptive_stats_t), intent(inout) :: stats
        real :: stddev_value

        if (stats%nr_values < 2) then
            write (unit=error_unit, fmt='(A)') &
                'error: too few data values to compute stddev'
            stop 2
        end if
        associate (n => stats%nr_values, sum => stats%sum, sum2 => stats%sum2)
            stddev_value = sqrt((sum2 - sum**2/n)/(n - 1.0))
        end associate
    end function get_stddev

    function get_median(stats) result(median_value)
        use, intrinsic :: iso_fortran_env, only : error_unit
        implicit none
        type(descriptive_stats_t), intent(inout) :: stats
        real :: median_value

        if (stats%nr_values == 0) then
            write (unit=error_unit, fmt='(A)') &
                'error: too few data value to compute median'
            stop 1
        end if
        associate(n => stats%nr_values, values => stats%values)
            call qsort(values)
            if (mod(n, 2) == 0) then
                median_value = (values(n/2) + values(n/2 + 1))/2.0
            else
                median_value = values(n/2 + 1)
            end if
        end associate
    end function get_median

    subroutine finalize_stats(stats)
        use, intrinsic :: iso_fortran_env, only : error_unit
        type(descriptive_stats_t), intent(inout) :: stats

        deallocate(stats%values)
    end subroutine finalize_stats

end module median_stats_mod

max_stars_repo_name: caguerra/Fortran-MOOC

max_issues_repo_name: caguerra/Fortran-MOOC

max_forks_repo_name: caguerra/Fortran-MOOC


text:
module bc_state_outlet_steady_1dchar
#include <messenger.h>
    use mod_kinds,              only: rk,ik
    use mod_constants,          only: ZERO, ONE, TWO, HALF, ME, CYLINDRICAL, NO_ID
    use mod_fluid,              only: gam
    use type_mesh,              only: mesh_t
    use type_bc_state,          only: bc_state_t
    use type_bc_patch,          only: bc_patch_t
    use type_chidg_worker,      only: chidg_worker_t
    use type_properties,        only: properties_t
    use type_element_info,      only: element_info_t, element_info
    use mod_chidg_mpi,          only: IRANK
    use mod_interpolate,        only: interpolate_face_autodiff
    use mpi_f08,                only: MPI_REAL8, MPI_SUM, MPI_AllReduce, mpi_comm, &
                                      MPI_INTEGER, MPI_BCast
    use DNAD_D
    implicit none

    !>  Name: Outlet - Average Pressure
    !!
    !!  Handle perturbations from average using one-dimensional characteristic analysis.
    !!
    !!  Options:
    !!      : Average Pressure
    !!
    !!  Behavior:
    !!
    !!  References:
    !!
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/23/2018
    !!
    !----------------------------------------------------------------------------------------
    type, public, extends(bc_state_t) :: outlet_steady_1dchar_t

    contains

        procedure   :: init                 ! Set-up bc state with options/name etc.
        procedure   :: init_bc_coupling     ! Implement specialized initialization procedure
        procedure   :: compute_bc_state     ! boundary condition function implementation
        procedure   :: compute_averages

    end type outlet_steady_1dchar_t
    !****************************************************************************************

contains

    !>
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/20/2018
    !!
    !--------------------------------------------------------------------------------
    subroutine init(self)
        class(outlet_steady_1dchar_t),   intent(inout) :: self

        ! Set name, family
        call self%set_name('Outlet - Steady 1D Characteristics')
        call self%set_family('Outlet')

        ! Add functions
        call self%bcproperties%add('Average Pressure','Required')

    end subroutine init
    !********************************************************************************

    !>  Initialize boundary group coupling.
    !!
    !!  Call global coupling routine to initialize implicit coupling between each
    !!  element with every other element on the boundary, a result of averaging
    !!  and Fourier transform operations.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/18/2018
    !!
    !--------------------------------------------------------------------------------
    subroutine init_bc_coupling(self,mesh,group_ID,bc_comm)
        class(outlet_steady_1dchar_t),  intent(inout)   :: self
        type(mesh_t),                   intent(inout)   :: mesh
        integer(ik),                    intent(in)      :: group_ID
        type(mpi_comm),                 intent(in)      :: bc_comm

        call self%init_bc_coupling_global(mesh,group_ID,bc_comm)

    end subroutine init_bc_coupling
    !********************************************************************************

    !>  Compute averaged quantities over the face.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   3/31/2017
    !!
    !!
    !--------------------------------------------------------------------------------
    subroutine compute_averages(self,worker,bc_COMM, v1_avg, v2_avg, v3_avg, vn_avg, density_avg, p_avg)
        class(outlet_steady_1dchar_t),   intent(inout)   :: self
        type(chidg_worker_t),               intent(inout)   :: worker
        type(mpi_comm),                     intent(in)      :: bc_COMM
        type(AD_D),                         intent(inout)   :: v1_avg
        type(AD_D),                         intent(inout)   :: v2_avg
        type(AD_D),                         intent(inout)   :: v3_avg
        type(AD_D),                         intent(inout)   :: vn_avg
        type(AD_D),                         intent(inout)   :: density_avg
        type(AD_D),                         intent(inout)   :: p_avg

        type(AD_D)          :: p_integral, v1_integral, v2_integral, v3_integral, &
                               vn_integral, density_integral, &
                               face_density, face_v1, face_v2, face_v3, face_vn, face_p
        type(element_info_t)    :: coupled_element

        type(AD_D), allocatable,    dimension(:)    ::  &
            density, mom1, mom2, mom3, energy, p,       &
            v1, v2, v3, vn

        real(rk),   allocatable,    dimension(:)    :: weights, areas, r

        integer(ik) :: ipatch, iface_bc, idomain_l, ielement_l, iface, ierr, itime, &
                       idensity, imom1, imom2, imom3, ienergy, group_ID, patch_ID, face_ID, &
                       icoupled, idomain_g_coupled, idomain_l_coupled, ielement_g_coupled,  &
                       ielement_l_coupled, iface_coupled, dof_start_coupled, coupled_iface
        real(rk)    :: face_area, total_area

        ! Zero integrated quantities
        total_area = ZERO

        ! Get location on domain
        idomain_l  = worker%element_info%idomain_l
        ielement_l = worker%element_info%ielement_l
        iface      = worker%iface

        ! Get location on bc_patch_group
        group_ID = worker%mesh%domain(idomain_l)%faces(ielement_l,iface)%group_ID
        patch_ID = worker%mesh%domain(idomain_l)%faces(ielement_l,iface)%patch_ID
        face_ID  = worker%mesh%domain(idomain_l)%faces(ielement_l,iface)%face_ID

        ! Loop through coupled faces and compute their contribution to the average pressure
        do icoupled = 1,worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%ncoupled_elements(face_ID)

            ! Get solution
            idensity = 1
            imom1    = 2
            imom2    = 3
            imom3    = 4
            ienergy  = 5
            itime    = 1

            ! Get face info from coupled element we want to interpolate from
            idomain_g_coupled  = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%idomain_g( icoupled)
            idomain_l_coupled  = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%idomain_l( icoupled)
            ielement_g_coupled = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%ielement_g(icoupled)
            ielement_l_coupled = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%ielement_l(icoupled)
            iface_coupled      = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%iface(     icoupled)
            dof_start_coupled  = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%dof_start( icoupled)

            coupled_element = element_info(idomain_g         = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%idomain_g(icoupled),          &
                                           idomain_l         = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%idomain_l(icoupled),          &
                                           ielement_g        = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%ielement_g(icoupled),         &
                                           ielement_l        = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%ielement_l(icoupled),         &
                                           iproc             = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%proc(icoupled),               &
                                           pelem_ID          = NO_ID,                                                                                &
                                           coordinate_system = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%coordinate_system(icoupled),  &
                                           eqn_ID            = NO_ID,                                                                                &
                                           nfields           = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%nfields(icoupled),            &
                                           ntime             = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%ntime(icoupled),              &
                                           nterms_s          = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%nterms_s(icoupled),           &
                                           nterms_c          = 0,                                                                                &
                                           dof_start         = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%dof_start(icoupled),          &
                                           dof_local_start   = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%dof_local_start(icoupled),    &
                                           xdof_start        = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%xdof_start(icoupled),         &
                                           xdof_local_start  = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%xdof_local_start(icoupled),   &
                                           recv_comm         = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%recv_comm(icoupled),          &
                                           recv_domain       = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%recv_domain(icoupled),        &
                                           recv_element      = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%recv_element(icoupled),       &
                                           recv_dof          = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%recv_dof(icoupled),           &
                                           recv_xdof         = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%recv_xdof(icoupled))

            coupled_iface = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%iface(icoupled)

            ! Interpolate coupled element solution on face of coupled element
            density = interpolate_face_autodiff(worker%mesh,worker%solverdata%q,coupled_element, worker%function_info, coupled_iface, idensity, itime, 'value', ME)
            mom1    = interpolate_face_autodiff(worker%mesh,worker%solverdata%q,coupled_element, worker%function_info, coupled_iface, imom1,    itime, 'value', ME)
            mom2    = interpolate_face_autodiff(worker%mesh,worker%solverdata%q,coupled_element, worker%function_info, coupled_iface, imom2,    itime, 'value', ME)
            mom3    = interpolate_face_autodiff(worker%mesh,worker%solverdata%q,coupled_element, worker%function_info, coupled_iface, imom3,    itime, 'value', ME)
            energy  = interpolate_face_autodiff(worker%mesh,worker%solverdata%q,coupled_element, worker%function_info, coupled_iface, ienergy,  itime, 'value', ME)

            ! TODO: fix for parallel!
            if (worker%mesh%domain(idomain_l_coupled)%elems(ielement_l_coupled)%coordinate_system == CYLINDRICAL) then
                mom2 = mom2 / worker%mesh%domain(idomain_l_coupled)%elems(ielement_l_coupled)%interp_coords_def(:,1)
            end if
            !r = worker%coordinate('1','boundary')
            !if (worker%coordinate_system() == 'Cylindrical') then
            !    mom2 = mom2 / r
            !end if

            ! Compute velocities and pressure
            v1 = mom1/density
            v2 = mom2/density
            v3 = mom3/density
            p = (gam-ONE)*(energy - HALF*(mom1*mom1 + mom2*mom2 + mom3*mom3)/density )

            ! TODO: fix for parallel!
            vn = v1*worker%mesh%domain(idomain_l_coupled)%faces(ielement_l_coupled,iface_coupled)%unorm_def(:,1) + &
                 v2*worker%mesh%domain(idomain_l_coupled)%faces(ielement_l_coupled,iface_coupled)%unorm_def(:,2) + &
                 v3*worker%mesh%domain(idomain_l_coupled)%faces(ielement_l_coupled,iface_coupled)%unorm_def(:,3)

            ! TODO: fix for parallel!
            ! Get weights + areas
            weights   = worker%mesh%domain(idomain_l)%faces(ielement_l,iface)%basis_s%weights_face(iface_coupled)
            areas     = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%data(icoupled)%areas
            face_area = worker%mesh%bc_patch_group(group_ID)%patch(patch_ID)%coupling(face_ID)%data(icoupled)%total_area

            ! Integrate and contribute to average
            face_v1      = sum(v1      * areas * weights)
            face_v2      = sum(v2      * areas * weights)
            face_v3      = sum(v3      * areas * weights)
            face_vn      = sum(vn      * areas * weights)
            face_density = sum(density * areas * weights)
            face_p       = sum(p       * areas * weights)

            if (allocated(v1_integral%xp_ad_)) then
                v1_integral = v1_integral + face_v1
            else
                v1_integral = face_v1
            end if

            if (allocated(v2_integral%xp_ad_)) then
                v2_integral = v2_integral + face_v2
            else
                v2_integral = face_v2
            end if

            if (allocated(v3_integral%xp_ad_)) then
                v3_integral = v3_integral + face_v3
            else
                v3_integral = face_v3
            end if

            if (allocated(vn_integral%xp_ad_)) then
                vn_integral = vn_integral + face_vn
            else
                vn_integral = face_vn
            end if

            if (allocated(p_integral%xp_ad_)) then
                p_integral = p_integral + face_p
            else
                p_integral = face_p
            end if

            if (allocated(density_integral%xp_ad_)) then
                density_integral = density_integral + face_density
            else
                density_integral = face_density
            end if

            total_area = total_area + face_area

        end do !icoupled

        ! Compute average pressure:
        !   area-weighted pressure integral over the total area
        v1_avg      = v1_integral      / total_area
        v2_avg      = v2_integral      / total_area
        v3_avg      = v3_integral      / total_area
        vn_avg      = vn_integral      / total_area
        density_avg = density_integral / total_area
        p_avg       = p_integral       / total_area

    end subroutine compute_averages
    !*******************************************************************************************

    !>  Compute routine for Pressure Outlet boundary condition state function.
    !!
    !!  @author Nathan A. steady_1dchar
    !!  @date   2/3/2016
    !!
    !!  @param[in]      worker  Interface for geometry, cache, integration, etc.
    !!  @param[inout]   prop    properties_t object containing equations and material_t objects
    !!
    !-------------------------------------------------------------------------------------------
    subroutine compute_bc_state(self,worker,prop,bc_COMM)
        class(outlet_steady_1dchar_t),    intent(inout)   :: self
        type(chidg_worker_t),               intent(inout)   :: worker
        class(properties_t),                intent(inout)   :: prop
        type(mpi_comm),                     intent(in)      :: bc_COMM

        ! Storage at quadrature nodes
        type(AD_D), allocatable, dimension(:)   ::                                      &
            density_m,  mom1_m,  mom2_m,  mom3_m,  energy_m,                            &
            density_bc, mom1_bc, mom2_bc, mom3_bc, energy_bc,                           &
            grad1_density_m, grad1_mom1_m, grad1_mom2_m, grad1_mom3_m, grad1_energy_m,  &
            grad2_density_m, grad2_mom1_m, grad2_mom2_m, grad2_mom3_m, grad2_energy_m,  &
            grad3_density_m, grad3_mom1_m, grad3_mom2_m, grad3_mom3_m, grad3_energy_m,  &
            v1_bc, v2_bc, v3_bc, p_bc,                                                  &
            v1_m,  v2_m,  v3_m,  vn_m, p_m,                                             &
            v1_t,  v2_t,  v3_t,                                                         &
            ddensity_c, dv1_c, dv2_c, dv3_c, dvn_c, dp_c,                               &
            c1, c2, c3, c4, c5, ddensity, dp, dv1, dv2, dv3, dvn,                       &
            dv1_avg, dv2_avg, dv3_avg, vn_avg_1, vn_avg_2, vn_avg_3

        type(AD_D)  :: p_avg, v1_avg, v2_avg, v3_avg, vn_avg, density_avg, M_avg, c_avg, c4_1d,    &
                       ddensity_avg, dvn_avg, dp_avg

        type(AD_D),     allocatable, dimension(:)   ::  p_user, r
        integer :: i

        ! Get back pressure from function.
        p_user = self%bcproperties%compute('Average Pressure',worker%time(),worker%coords())

        ! Interpolate interior solution to face quadrature nodes
        density_m = worker%get_field('Density'   , 'value', 'face interior')
        mom1_m    = worker%get_field('Momentum-1', 'value', 'face interior')
        mom2_m    = worker%get_field('Momentum-2', 'value', 'face interior')
        mom3_m    = worker%get_field('Momentum-3', 'value', 'face interior')
        energy_m  = worker%get_field('Energy'    , 'value', 'face interior')

        grad1_density_m = worker%get_field('Density'   , 'grad1', 'face interior')
        grad2_density_m = worker%get_field('Density'   , 'grad2', 'face interior')
        grad3_density_m = worker%get_field('Density'   , 'grad3', 'face interior')

        grad1_mom1_m    = worker%get_field('Momentum-1', 'grad1', 'face interior')
        grad2_mom1_m    = worker%get_field('Momentum-1', 'grad2', 'face interior')
        grad3_mom1_m    = worker%get_field('Momentum-1', 'grad3', 'face interior')

        grad1_mom2_m    = worker%get_field('Momentum-2', 'grad1', 'face interior')
        grad2_mom2_m    = worker%get_field('Momentum-2', 'grad2', 'face interior')
        grad3_mom2_m    = worker%get_field('Momentum-2', 'grad3', 'face interior')

        grad1_mom3_m    = worker%get_field('Momentum-3', 'grad1', 'face interior')
        grad2_mom3_m    = worker%get_field('Momentum-3', 'grad2', 'face interior')
        grad3_mom3_m    = worker%get_field('Momentum-3', 'grad3', 'face interior')

        grad1_energy_m  = worker%get_field('Energy'    , 'grad1', 'face interior')
        grad2_energy_m  = worker%get_field('Energy'    , 'grad2', 'face interior')
        grad3_energy_m  = worker%get_field('Energy'    , 'grad3', 'face interior')

        ! Store boundary gradient state. Grad(Q_bc). Do this here, before we
        ! compute any transformations for cylindrical.
        call worker%store_bc_state('Density'   , grad1_density_m, 'grad1')
        call worker%store_bc_state('Density'   , grad2_density_m, 'grad2')
        call worker%store_bc_state('Density'   , grad3_density_m, 'grad3')

        call worker%store_bc_state('Momentum-1', grad1_mom1_m,    'grad1')
        call worker%store_bc_state('Momentum-1', grad2_mom1_m,    'grad2')
        call worker%store_bc_state('Momentum-1', grad3_mom1_m,    'grad3')

        call worker%store_bc_state('Momentum-2', grad1_mom2_m,    'grad1')
        call worker%store_bc_state('Momentum-2', grad2_mom2_m,    'grad2')
        call worker%store_bc_state('Momentum-2', grad3_mom2_m,    'grad3')

        call worker%store_bc_state('Momentum-3', grad1_mom3_m,    'grad1')
        call worker%store_bc_state('Momentum-3', grad2_mom3_m,    'grad2')
        call worker%store_bc_state('Momentum-3', grad3_mom3_m,    'grad3')

        call worker%store_bc_state('Energy'    , grad1_energy_m,  'grad1')
        call worker%store_bc_state('Energy'    , grad2_energy_m,  'grad2')
        call worker%store_bc_state('Energy'    , grad3_energy_m,  'grad3')

        ! Account for cylindrical. Get tangential momentum from angular momentum.
        r = worker%coordinate('1','boundary')
        if (worker%coordinate_system() == 'Cylindrical') then
            mom2_m = mom2_m / r
            grad1_mom2_m = (grad1_mom2_m/r) - mom2_m/r
            grad2_mom2_m = (grad2_mom2_m/r)
            grad3_mom2_m = (grad3_mom2_m/r)
        end if

        ! Update average pressure
        call self%compute_averages(worker,bc_COMM,v1_avg, v2_avg, v3_avg, vn_avg, density_avg, p_avg)
        c_avg = sqrt(gam*p_avg/density_avg)

        ! Compute velocities
        v1_m = mom1_m/density_m
        v2_m = mom2_m/density_m
        v3_m = mom3_m/density_m
        vn_m = v1_m*worker%unit_normal(1) + &
               v2_m*worker%unit_normal(2) + &
               v3_m*worker%unit_normal(3)

        ! Compute tangential velocity part
        v1_t = v1_m - vn_m*worker%unit_normal(1)
        v2_t = v2_m - vn_m*worker%unit_normal(2)
        v3_t = v3_m - vn_m*worker%unit_normal(3)

        ! Compute pressure from extrapolated data
        p_m = worker%get_field('Pressure', 'value', 'face interior')

        ! Compute update for average quantities
        c4_1d        = TWO*(p_user(1) - p_avg)
        ddensity_avg =  c4_1d/(TWO*c_avg*c_avg)
        dvn_avg      = -c4_1d/(TWO*density_avg*c_avg)
        dp_avg       =  HALF*c4_1d

        ! Compute the average update along physical coordinates
        dv1_avg = dvn_avg*worker%unit_normal(1)
        dv2_avg = dvn_avg*worker%unit_normal(2)
        dv3_avg = dvn_avg*worker%unit_normal(3)

        ! Compute perturbation from avg
        dv1      = v1_m      - v1_avg
        dv2      = v2_m      - v2_avg
        dv3      = v3_m      - v3_avg
        dvn      = vn_m      - vn_avg
        ddensity = density_m - density_avg
        dp       = p_m       - p_avg

        ! Compute 1D characteristics
        allocate(c1(size(dp)), c2(size(dp)), c3(size(dp)), c4(size(dp)), c5(size(dp)))
        do i = 1,size(dp)
            c1(i) = -c_avg*c_avg*ddensity(i)   +  dp(i)
            c4(i) =  density_avg*c_avg*dvn(i)  +  dp(i)
            c5(i) = -density_avg*c_avg*dvn(i)  +  dp(i)
        end do

        ! Compute update from characteristics for perturbation quantities: No contrib from c5
        allocate(ddensity_c(size(dp)), dv1_c(size(dp)), dv2_c(size(dp)), dv3_c(size(dp)), dvn_c(size(dp)), dp_c(size(dp)))
        do i = 1,size(dp)
            ddensity_c(i) = -c1(i)/(c_avg*c_avg)  +  c4(i)/(TWO*c_avg*c_avg)
            dvn_c(i)      =  c4(i)/(TWO*density_avg*c_avg)
            dp_c(i)       =  c4(i)/TWO
        end do

        ! Compute velocity perturbation along physical coordinates
        dv1_c = dvn_c*worker%unit_normal(1)
        dv2_c = dvn_c*worker%unit_normal(2)
        dv3_c = dvn_c*worker%unit_normal(3)

        ! Compute contribution to each direction from average normal velocity
        vn_avg_1 = vn_avg*worker%unit_normal(1)
        vn_avg_2 = vn_avg*worker%unit_normal(2)
        vn_avg_3 = vn_avg*worker%unit_normal(3)

        ! Construct boundary state from average and perturbations
        density_bc = density_m
        v1_bc = density_m
        v2_bc = density_m
        v3_bc = density_m
        p_bc = density_m
        do i = 1,size(dp)
            density_bc(i) = density_avg  +  ddensity_avg  +  ddensity_c(i)
            v1_bc(i)      = vn_avg_1(i)  +  dv1_avg(i)    +  dv1_c(i)      +   v1_t(i)
            v2_bc(i)      = vn_avg_2(i)  +  dv2_avg(i)    +  dv2_c(i)      +   v2_t(i)
            v3_bc(i)      = vn_avg_3(i)  +  dv3_avg(i)    +  dv3_c(i)      +   v3_t(i)
            p_bc(i)       = p_avg        +  dp_avg        +  dp_c(i)
        end do

        ! Form conserved variables
        density_bc = density_bc
        mom1_bc    = density_bc*v1_bc
        mom2_bc    = density_bc*v2_bc
        mom3_bc    = density_bc*v3_bc
        energy_bc  = p_bc/(gam - ONE) + (density_bc*HALF)*(v1_bc*v1_bc + &
                                                           v2_bc*v2_bc + &
                                                           v3_bc*v3_bc)

        ! Account for cylindrical. Convert tangential momentum back to angular momentum.
        if (worker%coordinate_system() == 'Cylindrical') then
            mom2_bc = mom2_bc * r
        end if

        ! Store boundary condition state. Q_bc
        call worker%store_bc_state('Density'   , density_bc, 'value')
        call worker%store_bc_state('Momentum-1', mom1_bc,    'value')
        call worker%store_bc_state('Momentum-2', mom2_bc,    'value')
        call worker%store_bc_state('Momentum-3', mom3_bc,    'value')
        call worker%store_bc_state('Energy'    , energy_bc,  'value')

    end subroutine compute_bc_state
    !****************************************************************************************

end module bc_state_outlet_steady_1dchar

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG