MODULE caldyn_unstructured_mod
  USE ISO_C_BINDING
  USE OMP_LIB
  USE data_unstructured_mod
  IMPLICIT NONE
  SAVE

CONTAINS

#define INDICES1 ij,l,iq,iedge,edge,ivertex,vertex,ij_left,ij_right
#define INDICES2 ij_up,ij_down,itrisk,edge_trisk,kup,kdown
#define DECLARE_INDICES INTEGER INDICES1,INDICES2
#define PHI_BOT(ij) Phi_bot
#define PHI_BOT_VAR 0.

#define BINDC_(thename) BIND(C, name=#thename)
#define BINDC(thename) BINDC_(dynamico_ ## thename)

#define DBL REAL(C_DOUBLE)
#define DOUBLE1(m) DBL, DIMENSION(m)
#define DOUBLE2(m,n) DBL, DIMENSION(m,n)
#define DOUBLE3(m,n,p) DBL, DIMENSION(m,n,p)

#define FIELD_PS      DOUBLE1(primal_num)
#define FIELD_MASS    DOUBLE2(llm, primal_num)
#define FIELD_Z       DOUBLE2(llm, dual_num)
#define FIELD_U       DOUBLE2(llm, edge_num)
#define FIELD_UL      DOUBLE2(llm+1, edge_num)
#define FIELD_THETA   DOUBLE3(llm, primal_num, nqdyn)
#define FIELD_GEOPOT  DOUBLE2(llm+1, primal_num)

#define HASNAN(field) (ANY(.NOT.ABS(field)<1e20))

!----------------------------- Non-Hydrostatic -----------------------------

SUBROUTINE compute_NH_geopot(tau,dummy, m_ik, m_il, theta, W_il, Phi_il)
  FIELD_MASS   :: m_ik, theta, p_ik, A_ik, C_ik   ! IN*2,LOCAL*3
  FIELD_GEOPOT :: m_il, W_il, Phi_il, Phi_star_il, R_il, x_il, B_il, D_il  ! IN,INOUT*2, LOCAL*5 
  DBL :: tau, dummy, gamma, rho_ij, X_ij, Y_ij, wil, tau2_g, g2, gm2, ml_g2, c2_mik
  DECLARE_INDICES
  INTEGER :: iter
#include "../kernels_unst/compute_NH_geopot.k90"
END SUBROUTINE compute_NH_geopot

SUBROUTINE compute_caldyn_slow_NH(u,rhodz,Phi,W, hflux,du,dPhi,dW)
  FIELD_U      :: u, hflux, du   ! IN, OUT, OUT
  FIELD_MASS   :: rhodz, berni   ! IN, LOCAL
  FIELD_GEOPOT :: Phi,W,dPhi,dW, w_il, gradPhi2  ! IN,IN, OUT,OUT, LOCAL 
  FIELD_UL     :: DePhil, v_el, G_el, F_el ! LOCAL
  DECLARE_INDICES
  DBL :: W_el, W2_el, gPhi2, dP, divG, u2, uu
#include "../kernels_unst/caldyn_slow_NH.k90"
END SUBROUTINE compute_caldyn_slow_NH

SUBROUTINE compute_caldyn_solver(tau,rhodz,theta,pk,geopot,W,dPhi,dW,du)
  DBL, INTENT(IN) :: tau
  FIELD_MASS   :: rhodz,pk,berni,pres    ! IN, OUT, LOCAL
  FIELD_THETA  :: theta                  ! IN
  FIELD_GEOPOT :: geopot,W,dPhi,dW, m_il ! INOUT,INOUT, OUT,OUT, LOCAL 
  FIELD_U      :: du                     ! OUT
  DECLARE_INDICES
  DBL :: X_ij, rho_ij, T_ij, gamma, Cvd, vreff
#include "../kernels_unst/caldyn_solver.k90"
END SUBROUTINE compute_caldyn_solver

SUBROUTINE compute_caldyn_vert_NH(mass,geopot,W,wflux, du,dPhi,dW)
  FIELD_MASS   :: mass, eta_dot, wcov, W_etadot
  FIELD_GEOPOT :: geopot,W,wflux,dPhi,dW
  FIELD_U      :: du
  DECLARE_INDICES
  DBL :: w_ij, wflux_ij
#include "../kernels_unst/caldyn_vert_NH.k90"
END SUBROUTINE compute_caldyn_vert_NH

!----------------------------- Hydrostatic -------------------------------

SUBROUTINE compute_geopot(rhodz,theta,ps,pk,geopot)
  FIELD_MASS  :: rhodz,pk   ! IN, OUT
  FIELD_THETA :: theta      ! IN
  FIELD_GEOPOT :: geopot    ! IN(l=1)/OUT(l>1)
  FIELD_PS     :: ps        ! OUT
  DECLARE_INDICES
  DBL :: gdz, ke, uu, chi, gv, exner_ik, temp_ik, p_ik, qv, Rmix
#include "../kernels_unst/compute_geopot.k90"
END SUBROUTINE compute_geopot
!
SUBROUTINE compute_caldyn_slow_hydro(rhodz,theta,u, berni,hflux,du)
  FIELD_MASS  :: rhodz,berni! IN
  FIELD_THETA :: theta      ! IN
  FIELD_U     :: u,hflux,du ! IN, OUT, OUT
  DECLARE_INDICES
  LOGICAL, PARAMETER :: zero=.TRUE.
  DBL :: ke, uu
#include "../kernels_unst/caldyn_slow_hydro.k90"
END SUBROUTINE compute_caldyn_slow_hydro

!---------------------------------- Generic ------------------------------

SUBROUTINE caldyn_vert(convm,rhodz,theta,u, dmass_col,wflux,dtheta_rhodz,du, wwuu)
  FIELD_PS    :: dmass_col
  FIELD_MASS  :: convm, rhodz
  FIELD_U     :: u,du
  FIELD_THETA :: dtheta_rhodz, theta
  FIELD_GEOPOT :: wflux
  FIELD_UL     :: wwuu
  DECLARE_INDICES
  DBL :: dF_deta, dFu_deta
  wwuu=0.
  !$OMP BARRIER
#include "../kernels_unst/caldyn_wflux.k90"
#include "../kernels_unst/caldyn_dmass.k90"
#include "../kernels_unst/caldyn_vert.k90"
END SUBROUTINE caldyn_vert

SUBROUTINE compute_coriolis(hflux,theta,qu,Ftheta, convm,dtheta_rhodz,du)
  FIELD_U     :: hflux, Ftheta, qu, du
  FIELD_MASS  :: convm
  FIELD_THETA :: theta, dtheta_rhodz
  DECLARE_INDICES
  DBL :: divF, du_trisk
#include "../kernels_unst/coriolis.k90"
END SUBROUTINE

SUBROUTINE compute_theta(mass_col,rhodz,theta_rhodz, theta)
  FIELD_PS :: mass_col
  FIELD_MASS :: rhodz
  FIELD_THETA :: theta, theta_rhodz
  DECLARE_INDICES
  DBL :: m
#include "../kernels_unst/theta.k90"
END SUBROUTINE

SUBROUTINE compute_pvort_only(rhodz,u,qv,qu)
  FIELD_MASS :: rhodz
  FIELD_U    :: u,qu
  FIELD_Z    :: qv
  DECLARE_INDICES
  DBL :: etav, hv
#include "../kernels_unst/pvort_only.k90"
END SUBROUTINE compute_pvort_only

SUBROUTINE compute_caldyn_fast(tau, pk,berni,theta,geopot, du,u)
  DBL, INTENT(IN) :: tau
  FIELD_MASS   :: pk,berni  ! INOUT, OUT
  FIELD_THETA  :: theta     ! IN
  FIELD_GEOPOT :: geopot    ! IN
  FIELD_U      :: u,du      ! INOUT,INOUT
  DECLARE_INDICES
  DBL          :: due

#include "../kernels_unst/caldyn_fast.k90"

END SUBROUTINE compute_caldyn_fast

!----------------------------- Unused -----------------------------

SUBROUTINE gradient(b,grad) BINDC(gradient)
  DOUBLE2(llm,primal_num) :: b
  DOUBLE2(llm,edge_num)  :: grad
  DECLARE_INDICES
#include "../kernels_unst/gradient.k90"
END SUBROUTINE
!
SUBROUTINE div(u,divu) BINDC(div)
  DOUBLE2(llm,primal_num) :: divu
  DOUBLE2(llm,edge_num)  :: u
  DECLARE_INDICES
  DBL :: div_ij
  !$OMP PARALLEL NUM_THREADS(nb_threads)
#include "../kernels_unst/div.k90"
  !$OMP END PARALLEL
END SUBROUTINE

END MODULE caldyn_unstructured_mod