MODULE theta2theta_rhodz_mod

CONTAINS
  
  SUBROUTINE theta_rhodz2theta(f_ps,f_theta_rhodz,f_theta)
  USE icosa
  IMPLICIT NONE
    TYPE(t_field), POINTER :: f_ps(:)
    TYPE(t_field), POINTER :: f_theta_rhodz(:)
    TYPE(t_field), POINTER :: f_theta(:)
  
    REAL(rstd), POINTER :: ps(:)
    REAL(rstd), POINTER :: theta_rhodz(:,:)
    REAL(rstd), POINTER :: theta(:,:)
    INTEGER :: ind

    DO ind=1,ndomain
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind)
      ps=f_ps(ind)
      theta_rhodz=f_theta_rhodz(ind)
      theta=f_theta(ind)
      CALL compute_theta_rhodz2theta(ps, theta_rhodz,theta,0)
    ENDDO
  
  END SUBROUTINE theta_rhodz2theta

  SUBROUTINE theta_rhodz2temperature(f_ps,f_theta_rhodz,f_temp)
  USE icosa
  IMPLICIT NONE
    TYPE(t_field), POINTER :: f_ps(:)
    TYPE(t_field), POINTER :: f_theta_rhodz(:)
    TYPE(t_field), POINTER :: f_temp(:)
  
    REAL(rstd), POINTER :: ps(:)
    REAL(rstd), POINTER :: theta_rhodz(:,:)
    REAL(rstd), POINTER :: temp(:,:)
    INTEGER :: ind

    DO ind=1,ndomain
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind)
      ps=f_ps(ind)
      theta_rhodz=f_theta_rhodz(ind)
      temp=f_temp(ind)
      CALL compute_theta_rhodz2temperature(ps, theta_rhodz,temp,0)
    ENDDO
  
  END SUBROUTINE theta_rhodz2temperature
    
  SUBROUTINE theta2theta_rhodz(f_ps,f_theta,f_theta_rhodz)
  USE icosa
  IMPLICIT NONE
    TYPE(t_field), POINTER :: f_ps(:)
    TYPE(t_field), POINTER :: f_theta(:)
    TYPE(t_field), POINTER :: f_theta_rhodz(:)
  
    REAL(rstd), POINTER :: ps(:)
    REAL(rstd), POINTER :: theta(:,:)
    REAL(rstd), POINTER :: theta_rhodz(:,:)
    INTEGER :: ind

    DO ind=1,ndomain
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind)
      ps=f_ps(ind)
      theta=f_theta(ind)
      theta_rhodz=f_theta_rhodz(ind)
      CALL compute_theta2theta_rhodz(ps, theta, theta_rhodz,0)
    ENDDO
  
  END SUBROUTINE theta2theta_rhodz
  
  SUBROUTINE compute_theta2theta_rhodz(ps,theta, theta_rhodz,offset)
  USE icosa
  USE pression_mod
  IMPLICIT NONE
    REAL(rstd),INTENT(IN) :: ps(iim*jjm)
    REAL(rstd),INTENT(IN) :: theta(iim*jjm,llm)
    REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
    INTEGER,INTENT(IN) :: offset
    INTEGER :: i,j,ij,l
    REAL(rstd),ALLOCATABLE,SAVE :: p(:,:)

!$OMP BARRIER
!$OMP MASTER
      ALLOCATE( p(iim*jjm,llm+1))
!$OMP END MASTER
!$OMP BARRIER

    CALL compute_pression(ps,p,offset)
    
    DO    l    = 1, llm
!$OMP DO
      DO j=jj_begin-offset,jj_end+offset
        DO i=ii_begin-offset,ii_end+offset
          ij=(j-1)*iim+i
          theta_rhodz(ij,l) = theta(ij,l) * (p(ij,l)-p(ij,l+1))/g
        ENDDO
      ENDDO
    ENDDO

!$OMP BARRIER
!$OMP MASTER
      DEALLOCATE( p)
!$OMP END MASTER
!$OMP BARRIER    

  END SUBROUTINE compute_theta2theta_rhodz

  SUBROUTINE compute_theta_rhodz2theta(ps,theta_rhodz,theta,offset)
  USE icosa
  USE pression_mod
  IMPLICIT NONE
    REAL(rstd),INTENT(IN) :: ps(iim*jjm)
    REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm)
    REAL(rstd),INTENT(OUT) :: theta(iim*jjm,llm)
    INTEGER,INTENT(IN) :: offset
    INTEGER :: i,j,ij,l
    REAL(rstd),SAVE,ALLOCATABLE :: p(:,:)

!$OMP BARRIER
!$OMP MASTER
      ALLOCATE( p(iim*jjm,llm+1))
!$OMP END MASTER
!$OMP BARRIER
    CALL compute_pression(ps,p,offset)
    
    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
          theta(ij,l) = theta_rhodz(ij,l) / ((p(ij,l)-p(ij,l+1))/g)
        ENDDO
      ENDDO
    ENDDO

!$OMP BARRIER
!$OMP MASTER
      DEALLOCATE( p)
!$OMP END MASTER
!$OMP BARRIER    
    
  END SUBROUTINE compute_theta_rhodz2theta

  SUBROUTINE compute_theta_rhodz2temperature(ps,theta_rhodz,temp,offset)
  USE icosa
  USE pression_mod
  USE exner_mod
  IMPLICIT NONE
    REAL(rstd),INTENT(IN) :: ps(iim*jjm)
    REAL(rstd),INTENT(IN) :: theta_rhodz(iim*jjm,llm)
    REAL(rstd),INTENT(OUT) :: temp(iim*jjm,llm)
    INTEGER,INTENT(IN) :: offset
    INTEGER :: i,j,ij,l
    REAL(rstd) :: p(iim*jjm,llm+1)
    REAL(rstd) :: pk(iim*jjm,llm)
    REAL(rstd) :: pks(iim*jjm)

    CALL compute_pression(ps,p,offset)
    CALL compute_exner(ps,p,pks,pk,offset)
        
    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
          temp(ij,l) = ( theta_rhodz(ij,l) / ((p(ij,l)-p(ij,l+1))/g) ) * pk(ij,l) / cpp
        ENDDO
      ENDDO
    ENDDO
    
  END SUBROUTINE compute_theta_rhodz2temperature

  SUBROUTINE compute_temperature2theta_rhodz(ps,temp,theta_rhodz,offset)
  USE icosa
  USE pression_mod
  USE exner_mod
  IMPLICIT NONE
    REAL(rstd),INTENT(IN) :: ps(iim*jjm)
    REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
    REAL(rstd),INTENT(IN) :: temp(iim*jjm,llm)
    INTEGER,INTENT(IN) :: offset
    INTEGER :: i,j,ij,l
    REAL(rstd) :: p(iim*jjm,llm+1)
    REAL(rstd) :: pk(iim*jjm,llm)
    REAL(rstd) :: pks(iim*jjm)

    CALL compute_pression(ps,p,offset)
    CALL compute_exner(ps,p,pks,pk,offset)
        
    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
!          temp(ij,l) = ( theta_rhodz(ij,l) / ((p(ij,l)-p(ij,l+1))/g) ) * pk(ij,l) / cpp
          theta_rhodz(ij,l) = temp(ij,l) * ((p(ij,l)-p(ij,l+1))/g) / (pk(ij,l) / cpp )
        ENDDO
      ENDDO
    ENDDO
    
  END SUBROUTINE compute_temperature2theta_rhodz
   
END MODULE theta2theta_rhodz_mod