MODULE compute_caldyn_solver_mod
  USE grid_param, ONLY : llm
  IMPLICIT NONE
  PRIVATE

  PUBLIC :: compute_caldyn_solver

CONTAINS

  SUBROUTINE compute_caldyn_solver(tau,phis, rhodz,theta,pk, geopot,W, m_il,pres, dPhi,dW,du)
    USE icosa
    USE caldyn_vars_mod
    USE trace
    USE omp_para, ONLY : ll_begin, ll_end,ll_beginp1,ll_endp1
    USE disvert_mod, ONLY : ptop
    USE caldyn_kernels_hevi_mod
    REAL(rstd), PARAMETER :: pbot=1e5, rho_bot=1e6
    REAL(rstd),INTENT(IN) :: tau ! "solve" Phi-tau*dPhi/dt = Phi_rhs
    REAL(rstd),INTENT(IN)    :: phis(iim*jjm)
    REAL(rstd),INTENT(IN)    :: rhodz(iim*jjm,llm)
    REAL(rstd),INTENT(IN)    :: theta(iim*jjm,llm,nqdyn)
    REAL(rstd),INTENT(OUT)   :: pk(iim*jjm,llm)
    REAL(rstd),INTENT(INOUT) :: geopot(iim*jjm,llm+1)
    REAL(rstd),INTENT(INOUT) :: W(iim*jjm,llm+1) ! OUT if tau>0
    REAL(rstd),INTENT(OUT)   :: m_il(iim*jjm,llm+1)        ! rhodz averaged to interfaces
    REAL(rstd),INTENT(OUT)   :: pres(iim*jjm,llm)          ! pressure
    REAL(rstd),INTENT(OUT)   :: dW(iim*jjm,llm+1)
    REAL(rstd),INTENT(OUT)   :: dPhi(iim*jjm,llm+1)
    REAL(rstd),INTENT(OUT)   :: du(3*iim*jjm,llm)

    REAL(rstd) :: berni(iim*jjm,llm)         ! (W/m_il)^2
    REAL(rstd) :: berni1(iim*jjm)         ! (W/m_il)^2
    !REAL(rstd) :: gamma, rho_ij, T_ij, X_ij, Y_ij, vreff, Rd, Cvd, Rd_preff
    REAL(rstd) :: gamma, rho_ij, T_ij, X_ij, Y_ij, vreff, Cvd, Rd_preff
    INTEGER    :: ij, l

    CALL trace_start("compute_caldyn_solver")

    !Rd=cpp*kappa

    IF(dysl) THEN

!$OMP BARRIER

#include "../kernels_hex/caldyn_mil.k90"
  IF(tau>0) THEN ! solve implicit problem for geopotential
    CALL compute_NH_geopot(tau,phis, rhodz, m_il, theta, W, geopot)
  END IF
#define PHI_BOT(ij) phis(ij)
#include "../kernels_hex/caldyn_solver.k90"
#undef PHI_BOT
!$OMP BARRIER

    ELSE

#define BERNI(ij) berni1(ij)
    ! FIXME : vertical OpenMP parallelism will not work

    ! average m_ik to interfaces => m_il
    !DIR$ SIMD
    DO ij=ij_begin_ext,ij_end_ext
       m_il(ij,1) = .5*rhodz(ij,1)
    ENDDO
    DO l=2,llm
       !DIR$ SIMD
       DO ij=ij_begin_ext,ij_end_ext
          m_il(ij,l) = .5*(rhodz(ij,l-1)+rhodz(ij,l))
       ENDDO
    ENDDO
    !DIR$ SIMD
    DO ij=ij_begin_ext,ij_end_ext
       m_il(ij,llm+1) = .5*rhodz(ij,llm)
    ENDDO

    IF(tau>0) THEN ! solve implicit problem for geopotential
       CALL compute_NH_geopot(tau, phis, rhodz, m_il, theta, W, geopot)
    END IF

    ! Compute pressure, stored temporarily in pk
    ! kappa = R/Cp
    ! 1-kappa = Cv/Cp
    ! Cp/Cv = 1/(1-kappa)
    gamma = 1./(1.-kappa)
    DO l=1,llm
       !DIR$ SIMD
       DO ij=ij_begin_ext,ij_end_ext
          rho_ij = (g*rhodz(ij,l))/(geopot(ij,l+1)-geopot(ij,l))
          X_ij = (cpp/preff)*kappa*theta(ij,l,1)*rho_ij
          ! kappa.theta.rho = p/exner 
          ! => X = (p/p0)/(exner/Cp)
          !      = (p/p0)^(1-kappa)
          pk(ij,l) = preff*(X_ij**gamma)
       ENDDO
    ENDDO

    ! Update W, compute tendencies
    DO l=2,llm
       !DIR$ SIMD
       DO ij=ij_begin_ext,ij_end_ext
          dW(ij,l)   = (1./g)*(pk(ij,l-1)-pk(ij,l)) - m_il(ij,l)
          W(ij,l)    = W(ij,l)+tau*dW(ij,l) ! update W
          dPhi(ij,l) = g*g*W(ij,l)/m_il(ij,l)
       ENDDO
       !       PRINT *,'Max dPhi', l,ij_begin,ij_end, MAXVAL(abs(dPhi(ij_begin:ij_end,l)))
       !       PRINT *,'Max dW', l,ij_begin,ij_end, MAXVAL(abs(dW(ij_begin:ij_end,l))) 
    ENDDO
    ! Lower BC (FIXME : no orography yet !)
    DO ij=ij_begin,ij_end         
       dPhi(ij,1)=0
       W(ij,1)=0
       dW(ij,1)=0
       dPhi(ij,llm+1)=0 ! rigid lid
       W(ij,llm+1)=0
       dW(ij,llm+1)=0
    ENDDO
    ! Upper BC p=ptop
    !    DO ij=ij_omp_begin_ext,ij_omp_end_ext
    !       dPhi(ij,llm+1) = W(ij,llm+1)/rhodz(ij,llm)
    !       dW(ij,llm+1) = (1./g)*(pk(ij,llm)-ptop) - .5*rhodz(ij,llm)
    !    ENDDO
    
    ! Compute Exner function (needed by compute_caldyn_fast) and du=-g^2.grad(w^2)
    DO l=1,llm
       !DIR$ SIMD
       DO ij=ij_begin_ext,ij_end_ext
          pk(ij,l) = cpp*((pk(ij,l)/preff)**kappa) ! other formulae possible if exponentiation is slow
          BERNI(ij) = (-.25*g*g)*(              &
                 (W(ij,l)/m_il(ij,l))**2       &
               + (W(ij,l+1)/m_il(ij,l+1))**2 )
       ENDDO
       DO ij=ij_begin,ij_end         
          du(ij+u_right,l) = ne_right*(BERNI(ij)-BERNI(ij+t_right))
          du(ij+u_lup,l)   = ne_lup  *(BERNI(ij)-BERNI(ij+t_lup))
          du(ij+u_ldown,l) = ne_ldown*(BERNI(ij)-BERNI(ij+t_ldown))
       ENDDO
    ENDDO
#undef BERNI

    END IF ! dysl

    CALL trace_end("compute_caldyn_solver")
    
  END SUBROUTINE compute_caldyn_solver
  
END MODULE compute_caldyn_solver_mod