#include "cctk.h"

module qlm_variables
  use classify
  implicit none
  public
  save
  CCTK_REAL, allocatable, dimension(:,:) :: &
       qlm_gxx, qlm_gxy, qlm_gxz, qlm_gyy, qlm_gyz, qlm_gzz, &
       qlm_dgxxx, qlm_dgxyx, qlm_dgxzx, qlm_dgyyx, qlm_dgyzx, qlm_dgzzx, &
       qlm_dgxxy, qlm_dgxyy, qlm_dgxzy, qlm_dgyyy, qlm_dgyzy, qlm_dgzzy, &
       qlm_dgxxz, qlm_dgxyz, qlm_dgxzz, qlm_dgyyz, qlm_dgyzz, qlm_dgzzz, &
       qlm_ddgxxxx, qlm_ddgxyxx, qlm_ddgxzxx, qlm_ddgyyxx, qlm_ddgyzxx, qlm_ddgzzxx, &
       qlm_ddgxxxy, qlm_ddgxyxy, qlm_ddgxzxy, qlm_ddgyyxy, qlm_ddgyzxy, qlm_ddgzzxy, &
       qlm_ddgxxxz, qlm_ddgxyxz, qlm_ddgxzxz, qlm_ddgyyxz, qlm_ddgyzxz, qlm_ddgzzxz, &
       qlm_ddgxxyy, qlm_ddgxyyy, qlm_ddgxzyy, qlm_ddgyyyy, qlm_ddgyzyy, qlm_ddgzzyy, &
       qlm_ddgxxyz, qlm_ddgxyyz, qlm_ddgxzyz, qlm_ddgyyyz, qlm_ddgyzyz, qlm_ddgzzyz, &
       qlm_ddgxxzz, qlm_ddgxyzz, qlm_ddgxzzz, qlm_ddgyyzz, qlm_ddgyzzz, qlm_ddgzzzz, &
       qlm_kxx, qlm_kxy, qlm_kxz, qlm_kyy, qlm_kyz, qlm_kzz, &
       qlm_dkxxx, qlm_dkxyx, qlm_dkxzx, qlm_dkyyx, qlm_dkyzx, qlm_dkzzx, &
       qlm_dkxxy, qlm_dkxyy, qlm_dkxzy, qlm_dkyyy, qlm_dkyzy, qlm_dkzzy, &
       qlm_dkxxz, qlm_dkxyz, qlm_dkxzz, qlm_dkyyz, qlm_dkyzz, qlm_dkzzz, &
       qlm_alpha, &
       qlm_betax, qlm_betay, qlm_betaz, &
       qlm_ttt, &
       qlm_ttx, qlm_tty, qlm_ttz, &
       qlm_txx, qlm_txy, qlm_txz, qlm_tyy, qlm_tyz, qlm_tzz
  CCTK_REAL, allocatable, dimension(:,:) :: &
       qlm_dqttt, qlm_dqtpt, qlm_dqppt, qlm_dqttp, qlm_dqtpp, qlm_dqppp
  
  CCTK_REAL, allocatable, dimension(:,:) :: &
       qlm_spin_density
  
  type tetrad_derivs
     ! nabla_ll(a,b) = D_b l_a
     CCTK_REAL    :: nabla_ll(0:3,0:3)
     CCTK_REAL    :: nabla_nn(0:3,0:3)
     CCTK_COMPLEX :: nabla_mm(0:3,0:3)
  end type tetrad_derivs
  type(tetrad_derivs), allocatable, dimension(:,:) :: qlm_tetrad_derivs
  
contains
  
  subroutine allocate_variables(ntheta, nphi)
    integer,   intent(in) :: ntheta, nphi
    CCTK_REAL, parameter  :: zero = 0
    integer,   parameter  :: rk = kind(zero)
    type(tetrad_derivs)   :: tetrad_derivs_nan
    
    allocate(qlm_gxx(ntheta,nphi))
    allocate(qlm_gxy(ntheta,nphi))
    allocate(qlm_gxz(ntheta,nphi))
    allocate(qlm_gyy(ntheta,nphi))
    allocate(qlm_gyz(ntheta,nphi))
    allocate(qlm_gzz(ntheta,nphi))
    allocate(qlm_dgxxx(ntheta,nphi))
    allocate(qlm_dgxyx(ntheta,nphi))
    allocate(qlm_dgxzx(ntheta,nphi))
    allocate(qlm_dgyyx(ntheta,nphi))
    allocate(qlm_dgyzx(ntheta,nphi))
    allocate(qlm_dgzzx(ntheta,nphi))
    allocate(qlm_dgxxy(ntheta,nphi))
    allocate(qlm_dgxyy(ntheta,nphi))
    allocate(qlm_dgxzy(ntheta,nphi))
    allocate(qlm_dgyyy(ntheta,nphi))
    allocate(qlm_dgyzy(ntheta,nphi))
    allocate(qlm_dgzzy(ntheta,nphi))
    allocate(qlm_dgxxz(ntheta,nphi))
    allocate(qlm_dgxyz(ntheta,nphi))
    allocate(qlm_dgxzz(ntheta,nphi))
    allocate(qlm_dgyyz(ntheta,nphi))
    allocate(qlm_dgyzz(ntheta,nphi))
    allocate(qlm_dgzzz(ntheta,nphi))
    allocate(qlm_ddgxxxx(ntheta,nphi))
    allocate(qlm_ddgxyxx(ntheta,nphi))
    allocate(qlm_ddgxzxx(ntheta,nphi))
    allocate(qlm_ddgyyxx(ntheta,nphi))
    allocate(qlm_ddgyzxx(ntheta,nphi))
    allocate(qlm_ddgzzxx(ntheta,nphi))
    allocate(qlm_ddgxxxy(ntheta,nphi))
    allocate(qlm_ddgxyxy(ntheta,nphi))
    allocate(qlm_ddgxzxy(ntheta,nphi))
    allocate(qlm_ddgyyxy(ntheta,nphi))
    allocate(qlm_ddgyzxy(ntheta,nphi))
    allocate(qlm_ddgzzxy(ntheta,nphi))
    allocate(qlm_ddgxxxz(ntheta,nphi))
    allocate(qlm_ddgxyxz(ntheta,nphi))
    allocate(qlm_ddgxzxz(ntheta,nphi))
    allocate(qlm_ddgyyxz(ntheta,nphi))
    allocate(qlm_ddgyzxz(ntheta,nphi))
    allocate(qlm_ddgzzxz(ntheta,nphi))
    allocate(qlm_ddgxxyy(ntheta,nphi))
    allocate(qlm_ddgxyyy(ntheta,nphi))
    allocate(qlm_ddgxzyy(ntheta,nphi))
    allocate(qlm_ddgyyyy(ntheta,nphi))
    allocate(qlm_ddgyzyy(ntheta,nphi))
    allocate(qlm_ddgzzyy(ntheta,nphi))
    allocate(qlm_ddgxxyz(ntheta,nphi))
    allocate(qlm_ddgxyyz(ntheta,nphi))
    allocate(qlm_ddgxzyz(ntheta,nphi))
    allocate(qlm_ddgyyyz(ntheta,nphi))
    allocate(qlm_ddgyzyz(ntheta,nphi))
    allocate(qlm_ddgzzyz(ntheta,nphi))
    allocate(qlm_ddgxxzz(ntheta,nphi))
    allocate(qlm_ddgxyzz(ntheta,nphi))
    allocate(qlm_ddgxzzz(ntheta,nphi))
    allocate(qlm_ddgyyzz(ntheta,nphi))
    allocate(qlm_ddgyzzz(ntheta,nphi))
    allocate(qlm_ddgzzzz(ntheta,nphi))
    allocate(qlm_kxx(ntheta,nphi))
    allocate(qlm_kxy(ntheta,nphi))
    allocate(qlm_kxz(ntheta,nphi))
    allocate(qlm_kyy(ntheta,nphi))
    allocate(qlm_kyz(ntheta,nphi))
    allocate(qlm_kzz(ntheta,nphi))
    allocate(qlm_dkxxx(ntheta,nphi))
    allocate(qlm_dkxyx(ntheta,nphi))
    allocate(qlm_dkxzx(ntheta,nphi))
    allocate(qlm_dkyyx(ntheta,nphi))
    allocate(qlm_dkyzx(ntheta,nphi))
    allocate(qlm_dkzzx(ntheta,nphi))
    allocate(qlm_dkxxy(ntheta,nphi))
    allocate(qlm_dkxyy(ntheta,nphi))
    allocate(qlm_dkxzy(ntheta,nphi))
    allocate(qlm_dkyyy(ntheta,nphi))
    allocate(qlm_dkyzy(ntheta,nphi))
    allocate(qlm_dkzzy(ntheta,nphi))
    allocate(qlm_dkxxz(ntheta,nphi))
    allocate(qlm_dkxyz(ntheta,nphi))
    allocate(qlm_dkxzz(ntheta,nphi))
    allocate(qlm_dkyyz(ntheta,nphi))
    allocate(qlm_dkyzz(ntheta,nphi))
    allocate(qlm_dkzzz(ntheta,nphi))
    allocate(qlm_alpha(ntheta,nphi))
    allocate(qlm_betax(ntheta,nphi))
    allocate(qlm_betay(ntheta,nphi))
    allocate(qlm_betaz(ntheta,nphi))
    allocate(qlm_ttt(ntheta,nphi))
    allocate(qlm_ttx(ntheta,nphi))
    allocate(qlm_tty(ntheta,nphi))
    allocate(qlm_ttz(ntheta,nphi))
    allocate(qlm_txx(ntheta,nphi))
    allocate(qlm_txy(ntheta,nphi))
    allocate(qlm_txz(ntheta,nphi))
    allocate(qlm_tyy(ntheta,nphi))
    allocate(qlm_tyz(ntheta,nphi))
    allocate(qlm_tzz(ntheta,nphi))
    
    allocate(qlm_dqttt(ntheta,nphi))
    allocate(qlm_dqtpt(ntheta,nphi))
    allocate(qlm_dqppt(ntheta,nphi))
    allocate(qlm_dqttp(ntheta,nphi))
    allocate(qlm_dqtpp(ntheta,nphi))
    allocate(qlm_dqppp(ntheta,nphi))
    
    allocate(qlm_tetrad_derivs(ntheta,nphi))
    
    allocate(qlm_spin_density(ntheta,nphi))
    
    qlm_gxx = TAT_nan()
    qlm_gxy = TAT_nan()
    qlm_gxz = TAT_nan()
    qlm_gyy = TAT_nan()
    qlm_gyz = TAT_nan()
    qlm_gzz = TAT_nan()
    qlm_dgxxx = TAT_nan()
    qlm_dgxyx = TAT_nan()
    qlm_dgxzx = TAT_nan()
    qlm_dgyyx = TAT_nan()
    qlm_dgyzx = TAT_nan()
    qlm_dgzzx = TAT_nan()
    qlm_dgxxy = TAT_nan()
    qlm_dgxyy = TAT_nan()
    qlm_dgxzy = TAT_nan()
    qlm_dgyyy = TAT_nan()
    qlm_dgyzy = TAT_nan()
    qlm_dgzzy = TAT_nan()
    qlm_dgxxz = TAT_nan()
    qlm_dgxyz = TAT_nan()
    qlm_dgxzz = TAT_nan()
    qlm_dgyyz = TAT_nan()
    qlm_dgyzz = TAT_nan()
    qlm_dgzzz = TAT_nan()
    qlm_ddgxxxx = TAT_nan()
    qlm_ddgxyxx = TAT_nan()
    qlm_ddgxzxx = TAT_nan()
    qlm_ddgyyxx = TAT_nan()
    qlm_ddgyzxx = TAT_nan()
    qlm_ddgzzxx = TAT_nan()
    qlm_ddgxxxy = TAT_nan()
    qlm_ddgxyxy = TAT_nan()
    qlm_ddgxzxy = TAT_nan()
    qlm_ddgyyxy = TAT_nan()
    qlm_ddgyzxy = TAT_nan()
    qlm_ddgzzxy = TAT_nan()
    qlm_ddgxxxz = TAT_nan()
    qlm_ddgxyxz = TAT_nan()
    qlm_ddgxzxz = TAT_nan()
    qlm_ddgyyxz = TAT_nan()
    qlm_ddgyzxz = TAT_nan()
    qlm_ddgzzxz = TAT_nan()
    qlm_ddgxxyy = TAT_nan()
    qlm_ddgxyyy = TAT_nan()
    qlm_ddgxzyy = TAT_nan()
    qlm_ddgyyyy = TAT_nan()
    qlm_ddgyzyy = TAT_nan()
    qlm_ddgzzyy = TAT_nan()
    qlm_ddgxxyz = TAT_nan()
    qlm_ddgxyyz = TAT_nan()
    qlm_ddgxzyz = TAT_nan()
    qlm_ddgyyyz = TAT_nan()
    qlm_ddgyzyz = TAT_nan()
    qlm_ddgzzyz = TAT_nan()
    qlm_ddgxxzz = TAT_nan()
    qlm_ddgxyzz = TAT_nan()
    qlm_ddgxzzz = TAT_nan()
    qlm_ddgyyzz = TAT_nan()
    qlm_ddgyzzz = TAT_nan()
    qlm_ddgzzzz = TAT_nan()
    qlm_kxx = TAT_nan()
    qlm_kxy = TAT_nan()
    qlm_kxz = TAT_nan()
    qlm_kyy = TAT_nan()
    qlm_kyz = TAT_nan()
    qlm_kzz = TAT_nan()
    qlm_dkxxx = TAT_nan()
    qlm_dkxyx = TAT_nan()
    qlm_dkxzx = TAT_nan()
    qlm_dkyyx = TAT_nan()
    qlm_dkyzx = TAT_nan()
    qlm_dkzzx = TAT_nan()
    qlm_dkxxy = TAT_nan()
    qlm_dkxyy = TAT_nan()
    qlm_dkxzy = TAT_nan()
    qlm_dkyyy = TAT_nan()
    qlm_dkyzy = TAT_nan()
    qlm_dkzzy = TAT_nan()
    qlm_dkxxz = TAT_nan()
    qlm_dkxyz = TAT_nan()
    qlm_dkxzz = TAT_nan()
    qlm_dkyyz = TAT_nan()
    qlm_dkyzz = TAT_nan()
    qlm_dkzzz = TAT_nan()
    qlm_alpha = TAT_nan()
    qlm_betax = TAT_nan()
    qlm_betay = TAT_nan()
    qlm_betaz = TAT_nan()
    qlm_ttt = TAT_nan()
    qlm_ttx = TAT_nan()
    qlm_tty = TAT_nan()
    qlm_ttz = TAT_nan()
    qlm_txx = TAT_nan()
    qlm_txy = TAT_nan()
    qlm_txz = TAT_nan()
    qlm_tyy = TAT_nan()
    qlm_tyz = TAT_nan()
    qlm_tzz = TAT_nan()
    
    qlm_dqttt = TAT_nan()
    qlm_dqtpt = TAT_nan()
    qlm_dqppt = TAT_nan()
    qlm_dqttp = TAT_nan()
    qlm_dqtpp = TAT_nan()
    qlm_dqppp = TAT_nan()
    
    tetrad_derivs_nan%nabla_ll = TAT_nan()
    tetrad_derivs_nan%nabla_nn = TAT_nan()
    tetrad_derivs_nan%nabla_mm = cmplx(TAT_nan(),TAT_nan(),rk)
    qlm_tetrad_derivs = tetrad_derivs_nan

    qlm_spin_density = TAT_nan()
  end subroutine allocate_variables

  subroutine deallocate_variables
    deallocate(qlm_gxx)
    deallocate(qlm_gxy)
    deallocate(qlm_gxz)
    deallocate(qlm_gyy)
    deallocate(qlm_gyz)
    deallocate(qlm_gzz)
    deallocate(qlm_dgxxx)
    deallocate(qlm_dgxyx)
    deallocate(qlm_dgxzx)
    deallocate(qlm_dgyyx)
    deallocate(qlm_dgyzx)
    deallocate(qlm_dgzzx)
    deallocate(qlm_dgxxy)
    deallocate(qlm_dgxyy)
    deallocate(qlm_dgxzy)
    deallocate(qlm_dgyyy)
    deallocate(qlm_dgyzy)
    deallocate(qlm_dgzzy)
    deallocate(qlm_dgxxz)
    deallocate(qlm_dgxyz)
    deallocate(qlm_dgxzz)
    deallocate(qlm_dgyyz)
    deallocate(qlm_dgyzz)
    deallocate(qlm_dgzzz)
    deallocate(qlm_ddgxxxx)
    deallocate(qlm_ddgxyxx)
    deallocate(qlm_ddgxzxx)
    deallocate(qlm_ddgyyxx)
    deallocate(qlm_ddgyzxx)
    deallocate(qlm_ddgzzxx)
    deallocate(qlm_ddgxxxy)
    deallocate(qlm_ddgxyxy)
    deallocate(qlm_ddgxzxy)
    deallocate(qlm_ddgyyxy)
    deallocate(qlm_ddgyzxy)
    deallocate(qlm_ddgzzxy)
    deallocate(qlm_ddgxxxz)
    deallocate(qlm_ddgxyxz)
    deallocate(qlm_ddgxzxz)
    deallocate(qlm_ddgyyxz)
    deallocate(qlm_ddgyzxz)
    deallocate(qlm_ddgzzxz)
    deallocate(qlm_ddgxxyy)
    deallocate(qlm_ddgxyyy)
    deallocate(qlm_ddgxzyy)
    deallocate(qlm_ddgyyyy)
    deallocate(qlm_ddgyzyy)
    deallocate(qlm_ddgzzyy)
    deallocate(qlm_ddgxxyz)
    deallocate(qlm_ddgxyyz)
    deallocate(qlm_ddgxzyz)
    deallocate(qlm_ddgyyyz)
    deallocate(qlm_ddgyzyz)
    deallocate(qlm_ddgzzyz)
    deallocate(qlm_ddgxxzz)
    deallocate(qlm_ddgxyzz)
    deallocate(qlm_ddgxzzz)
    deallocate(qlm_ddgyyzz)
    deallocate(qlm_ddgyzzz)
    deallocate(qlm_ddgzzzz)
    deallocate(qlm_kxx)
    deallocate(qlm_kxy)
    deallocate(qlm_kxz)
    deallocate(qlm_kyy)
    deallocate(qlm_kyz)
    deallocate(qlm_kzz)
    deallocate(qlm_dkxxx)
    deallocate(qlm_dkxyx)
    deallocate(qlm_dkxzx)
    deallocate(qlm_dkyyx)
    deallocate(qlm_dkyzx)
    deallocate(qlm_dkzzx)
    deallocate(qlm_dkxxy)
    deallocate(qlm_dkxyy)
    deallocate(qlm_dkxzy)
    deallocate(qlm_dkyyy)
    deallocate(qlm_dkyzy)
    deallocate(qlm_dkzzy)
    deallocate(qlm_dkxxz)
    deallocate(qlm_dkxyz)
    deallocate(qlm_dkxzz)
    deallocate(qlm_dkyyz)
    deallocate(qlm_dkyzz)
    deallocate(qlm_dkzzz)
    deallocate(qlm_alpha)
    deallocate(qlm_betax)
    deallocate(qlm_betay)
    deallocate(qlm_betaz)
    deallocate(qlm_ttt)
    deallocate(qlm_ttx)
    deallocate(qlm_tty)
    deallocate(qlm_ttz)
    deallocate(qlm_txx)
    deallocate(qlm_txy)
    deallocate(qlm_txz)
    deallocate(qlm_tyy)
    deallocate(qlm_tyz)
    deallocate(qlm_tzz)
    
    deallocate(qlm_dqttt)
    deallocate(qlm_dqtpt)
    deallocate(qlm_dqppt)
    deallocate(qlm_dqttp)
    deallocate(qlm_dqtpp)
    deallocate(qlm_dqppp)
    
    deallocate(qlm_tetrad_derivs)

    deallocate(qlm_spin_density)
  end subroutine deallocate_variables

end module qlm_variables