MODULE observable_mod
  USE icosa
  IMPLICIT NONE
  PRIVATE

  TYPE(t_field),POINTER, SAVE :: f_buf_i(:), f_buf_ulon(:), f_buf_ulat(:), f_buf_u3d(:)
  TYPE(t_field),POINTER, SAVE :: f_buf1_i(:), f_buf2_i(:)
  TYPE(t_field),POINTER, SAVE :: f_buf_v(:), f_buf_s(:), f_buf_p(:)

! temporary shared variable for caldyn
  TYPE(t_field),POINTER, SAVE :: f_theta(:)

  PUBLIC init_observable, write_output_fields_basic, f_theta

CONTAINS
  
  SUBROUTINE init_observable
    CALL allocate_field(f_buf_i,   field_t,type_real,llm,name="buffer_i")
    CALL allocate_field(f_buf_p,   field_t,type_real,llm+1) 
    CALL allocate_field(f_buf_u3d, field_t,type_real,3,llm)  ! 3D vel at cell centers
    CALL allocate_field(f_buf_ulon,field_t,type_real,llm, name="buf_ulon")
    CALL allocate_field(f_buf_ulat,field_t,type_real,llm, name="buf_ulat")
    CALL allocate_field(f_buf_v,   field_z,type_real,llm)
    CALL allocate_field(f_buf_s,   field_t,type_real)

    CALL allocate_field(f_theta, field_t,type_real,llm,  name='theta')   ! potential temperature
  END SUBROUTINE init_observable
  
  SUBROUTINE write_output_fields_basic(f_ps, f_u, f_q)
    USE wind_mod
    USE output_field_mod
    USE omp_para
    TYPE(t_field),POINTER :: f_ps(:), f_u(:), f_q(:)
!    IF (is_master) PRINT *,'CALL write_output_fields_basic'
    CALL un2ulonlat(f_u, f_buf_ulon, f_buf_ulat)
    CALL output_field("ulon",f_buf_ulon)
    CALL output_field("ulat",f_buf_ulat)
    CALL output_field("ps",f_ps)
    CALL output_field("Ai",geom%Ai)
    !       CALL output_field("dps",f_dps)
    !CALL output_field("mass",f_mass)
    !       CALL output_field("dmass",f_dmass)
    !       CALL output_field("vort",f_qv)
    CALL output_field("theta",f_theta)
    !       CALL output_field("exner",f_pk)
    !       CALL output_field("pv",f_qv)
    CALL output_field("q",f_q)
  END SUBROUTINE write_output_fields_basic
  
  SUBROUTINE write_output_fields(f_ps, f_phis, f_dps, f_u, f_theta_rhodz, f_q, &
       f_buf_i, f_buf_v, f_buf_i3, f_buf1_i, f_buf2_i, f_buf_s, f_buf_p)
    USE icosa
    USE vorticity_mod
    USE theta2theta_rhodz_mod
    USE pression_mod
    USE omega_mod
    USE write_field_mod
    USE vertical_interp_mod
    USE wind_mod
    TYPE(t_field),POINTER :: f_ps(:), f_phis(:), f_u(:), f_theta_rhodz(:), f_q(:), f_dps(:), &
         f_buf_i(:), f_buf_v(:), f_buf_i3(:), f_buf1_i(:), f_buf2_i(:), f_buf_s(:), f_buf_p(:)
    
    REAL(rstd) :: out_pression_level
    CHARACTER(LEN=255) :: str_pression
    CHARACTER(LEN=255) :: physics_type
    
    out_pression_level=0.
    CALL getin("out_pression_level",out_pression_level) 
    WRITE(str_pression,*) INT(out_pression_level/100)
    str_pression=ADJUSTL(str_pression)
    
    CALL writefield("ps",f_ps)
    CALL writefield("dps",f_dps)
    CALL writefield("phis",f_phis)
    CALL vorticity(f_u,f_buf_v)
    CALL writefield("vort",f_buf_v)

    CALL w_omega(f_ps, f_u, f_buf_i)
    CALL writefield('omega', f_buf_i)
    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
      CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
      CALL writefield("omega"//TRIM(str_pression),f_buf_s)
    ENDIF
    
    ! Temperature
!    CALL theta_rhodz2temperature(f_ps,f_theta_rhodz,f_buf_i) ; ! FIXME
     
    CALL getin('physics',physics_type)
    IF (TRIM(physics_type)=='dcmip') THEN
      CALL Tv2T(f_buf_i,f_q,f_buf1_i) 
      CALL writefield("T",f_buf1_i)
      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
        CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
        CALL writefield("T"//TRIM(str_pression),f_buf_s)
      ENDIF
    ELSE
      CALL writefield("T",f_buf_i)
      IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
        CALL vertical_interp(f_ps,f_buf_i,f_buf_s,out_pression_level)
        CALL writefield("T"//TRIM(str_pression),f_buf_s)
      ENDIF
    ENDIF
   
    ! velocity components
    CALL un2ulonlat(f_u, f_buf1_i, f_buf2_i)
    CALL writefield("ulon",f_buf1_i)
    CALL writefield("ulat",f_buf2_i)

    IF (out_pression_level<=preff .AND. out_pression_level > 0) THEN
      CALL vertical_interp(f_ps,f_buf1_i,f_buf_s,out_pression_level)
      CALL writefield("ulon"//TRIM(str_pression),f_buf_s)
      CALL vertical_interp(f_ps,f_buf2_i,f_buf_s,out_pression_level)
      CALL writefield("ulat"//TRIM(str_pression),f_buf_s)
    ENDIF
    
    ! geopotential ! FIXME
    CALL thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_buf_s,f_buf_p,f_buf1_i,f_buf2_i,f_buf_i)
    CALL writefield("p",f_buf_p)
!    CALL writefield("phi",f_geopot)   ! geopotential
    CALL writefield("theta",f_buf1_i) ! potential temperature
    CALL writefield("pk",f_buf2_i)    ! Exner pressure
  
  END SUBROUTINE write_output_fields
  
  SUBROUTINE thetarhodz2geopot(f_ps,f_phis,f_theta_rhodz, f_pks,f_p,f_theta,f_pk,f_phi) 
    USE field_mod
    USE pression_mod
    USE exner_mod
    USE geopotential_mod
    USE theta2theta_rhodz_mod
    TYPE(t_field), POINTER :: f_ps(:), f_phis(:), f_theta_rhodz(:), &  ! IN
         f_pks(:), f_p(:), f_theta(:), f_pk(:), f_phi(:)               ! OUT
    REAL(rstd),POINTER :: pk(:,:), p(:,:), theta(:,:), theta_rhodz(:,:), &
         phi(:,:), phis(:), ps(:), pks(:)
    INTEGER :: ind
    
    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)
!$OMP BARRIER
       CALL compute_pression(ps,p,0)
       pk = f_pk(ind)
       pks = f_pks(ind)
!$OMP BARRIER
       CALL compute_exner(ps,p,pks,pk,0)
!$OMP BARRIER
       theta_rhodz = f_theta_rhodz(ind)
       theta = f_theta(ind)
       CALL compute_theta_rhodz2theta(ps, theta_rhodz,theta,0)
       phis = f_phis(ind)
       phi = f_phi(ind)
       CALL compute_geopotential(phis,pks,pk,theta,phi,0)
    END DO

  END SUBROUTINE thetarhodz2geopot
  
  SUBROUTINE Tv2T(f_Tv, f_q, f_T)
    TYPE(t_field), POINTER :: f_TV(:)
    TYPE(t_field), POINTER :: f_q(:)
    TYPE(t_field), POINTER :: f_T(:)
    
    REAL(rstd),POINTER :: Tv(:,:), q(:,:,:), T(:,:)
    INTEGER :: ind
    
    DO ind=1,ndomain
       IF (.NOT. assigned_domain(ind)) CYCLE
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
       Tv=f_Tv(ind)
       q=f_q(ind)
       T=f_T(ind)
       T=Tv/(1+0.608*q(:,:,1))
    END DO
    
  END SUBROUTINE Tv2T
  
END MODULE observable_mod