MODULE exner_mod

  INTEGER,SAVE :: caldyn_exner
!$OMP THREADPRIVATE(caldyn_exner)

  INTEGER, PARAMETER :: lmdz=3, direct=4

CONTAINS
  
  SUBROUTINE exner(f_ps,f_p,f_pks,f_pk)
  USE icosa
  IMPLICIT NONE
    TYPE(t_field), POINTER :: f_ps(:)  ! IN
    TYPE(t_field), POINTER :: f_p(:)   ! IN
    TYPE(t_field), POINTER :: f_pks(:) ! OUT
    TYPE(t_field), POINTER :: f_pk(:)  ! OUT
  
    REAL(rstd), POINTER :: ps(:)
    REAL(rstd), POINTER :: p(:,:)
    REAL(rstd), POINTER :: pks(:)
    REAL(rstd), POINTER :: pk(:,:)
    INTEGER :: ind

!$OMP BARRIER
    DO ind=1,ndomain
      IF (.NOT. assigned_domain(ind)) CYCLE
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind)
      ps=f_ps(ind)
      p=f_p(ind)
      pks=f_pks(ind)
      pk=f_pk(ind)
      CALL compute_exner(ps, p, pks, pk, 0)
    ENDDO
!$OMP BARRIER
  
  END SUBROUTINE exner
  
  SUBROUTINE compute_exner(ps,p,pks,pk,offset)
  USE icosa
  USE disvert_mod
  USE pression_mod
  USE omp_para
  IMPLICIT NONE
    REAL(rstd),INTENT(IN) :: ps(iim*jjm)
    REAL(rstd),INTENT(IN) :: p(iim*jjm,llm+1)
    REAL(rstd),INTENT(OUT) :: pks(iim*jjm)
    REAL(rstd),INTENT(OUT) :: pk(iim*jjm,llm)
    INTEGER,INTENT(IN) :: offset
    INTEGER :: i,j,ij,l
    
    ! surface : pks
    IF (is_omp_level_master) THEN
       
       DO j=jj_begin-offset,jj_end+offset
          DO i=ii_begin-offset,ii_end+offset
             ij=(j-1)*iim+i
             pks(ij) = cpp * ( ps(ij)/preff ) ** kappa
          ENDDO
       ENDDO
       
    ENDIF
    
    ! 3D : pk
    DO l = 1, llm
       DO j=jj_begin-offset,jj_end+offset
          DO i=ii_begin-offset,ii_end+offset
             ij=(j-1)*iim+i
             pk(ij,l) = cpp * ((.5/preff)*(p(ij,l)+p(ij,l+1))) ** kappa
             ENDDO
          ENDDO
       ENDDO
       
  END SUBROUTINE compute_exner
  
END MODULE exner_mod