module pcm_thermodynamics_mod

      use comcstfi_mod, only: cppd_ref,cppv_ref,rd_ref, rcp_ref, mugaz_ref
      
      implicit none
      
      REAL, SAVE, ALLOCATABLE, PROTECTED :: r(:,:)     ! Reduced gas constant of atmosphere
      REAL, SAVE, ALLOCATABLE, PROTECTED :: cpp(:,:)   ! Cp of the atmosphere (J/K/kg)
      REAL, SAVE, ALLOCATABLE, PROTECTED :: rcp(:,:)   ! r/cpp (gamma) of atmosphere
!$OMP THREADPRIVATE(r,cpp,rcp)
      
      contains
      
      subroutine thermodynamics_update(thermo_phy, pplay, pplev, t, ngrid, nlayer, nq, q, iq, zh, zpopsk)
        CHARACTER(LEN=*), INTENT(IN) :: thermo_phy               ! flag
        REAL, INTENT(IN)    :: pplay(ngrid,nlayer)    ! pressure (Pa) at mid-layer
        REAL, INTENT(IN)    :: pplev(ngrid,nlayer+1)  ! pressure (Pa) at layer interfaces
        REAL, INTENT(IN)    :: t(ngrid,nlayer)      ! temperature (K)
        INTEGER, INTENT(IN) :: ngrid, nlayer,nq     ! Number of cells, vertical layers and tracers
        REAL, INTENT(IN)    :: q(ngrid,nlayer,nq)   ! Mass mixing ratio of tracers (kg/kg_of_air)
        INTEGER, INTENT(IN) :: iq                   ! id variable active tracer
        REAL, INTENT(OUT)   :: zh(ngrid,nlayer)     ! Potential temperature (K)
        REAL, INTENT(OUT)   :: zpopsk(ngrid,nlayer) ! Exner function
    
        ! Local variables
        integer :: ig,l
        logical, save :: firstcall=.true.
!$OMP THREADPRIVATE(firstcall)
        character(len=80),parameter :: myname = "thermodynamics_update" 
        
        if (firstcall) then
          ALLOCATE(r(ngrid,nlayer))
          ALLOCATE(cpp(ngrid,nlayer))
          ALLOCATE(rcp(ngrid,nlayer))
          
          SELECT CASE (TRIM(thermo_phy))
            CASE('thermo_uni_ideal')
              ! Ideal gas, homogeneous
              r(:,:) = rd_ref
              cpp(:,:) = cppd_ref
              rcp(:,:) = rcp_ref
            CASE DEFAULT
            write(*,*) 'Bad selector for thermodynamics mod: <', TRIM(thermo_phy), '>'
            call abort_physic(trim(myname),'Bad selector for thermodynamics mod!',1)
          END SELECT
          
          firstcall=.false.
        endif ! of if (firstcall)
        
        SELECT CASE (TRIM(thermo_phy))

         CASE('thermo_uni_ideal')
            ! Ideal gas
            do l=1,nlayer
              do ig=1,ngrid
                zpopsk(ig,l) = (pplay(ig,l)/pplev(ig,1))**rcp(ig,l)
                zh(ig,l) = t(ig,l)/zpopsk(ig,l)
              enddo
            enddo
         CASE DEFAULT
            write(*,*) 'Bad selector for thermodynamics mod: <', TRIM(thermo_phy), '>'
            write(*,*) 'Option is <thermo_uni_ideal>'
            call abort_physic(trim(myname),'Bad selector for thermodynamics mod!',1)
        END SELECT
        
        end subroutine thermodynamics_update
        
        end module pcm_thermodynamics_mod