MODULE physics_external_mod
  USE field_mod
  
  INTEGER,SAVE :: it
!$OMP THREADPRIVATE(it)

  TYPE(t_field),POINTER,SAVE :: f_phis(:)
  TYPE(t_field),POINTER,SAVE :: f_ps(:)
  TYPE(t_field),POINTER,SAVE :: f_theta_rhodz(:)
  TYPE(t_field),POINTER,SAVE :: f_u(:)
  TYPE(t_field),POINTER,SAVE :: f_wflux(:)
  TYPE(t_field),POINTER,SAVE :: f_q(:)

  TYPE(t_field),POINTER,SAVE :: f_theta_rhodz0(:)
  TYPE(t_field),POINTER,SAVE :: f_u0(:)
  TYPE(t_field),POINTER,SAVE :: f_q0(:)
 
  TYPE(t_field),POINTER,SAVE :: f_dtheta_rhodz(:)
  TYPE(t_field),POINTER,SAVE :: f_du(:)
  TYPE(t_field),POINTER,SAVE :: f_dq(:)

  TYPE(t_field),POINTER,SAVE :: f_rhodz(:)
  TYPE(t_field),POINTER,SAVE :: f_rhodz0(:)
  
  LOGICAL,SAVE :: phys_smooth_tendency
!$OMP THREADPRIVATE(phys_smooth_tendency)  


CONTAINS

  SUBROUTINE init_physics
  USE icosa
  IMPLICIT NONE

    CALL initialize_external_physics
!$OMP PARALLEL    
    CALL allocate_field(f_theta_rhodz0, field_t, type_real, llm, nqdyn, name='theta_rhodz0')
    CALL allocate_field(f_u0,field_u,type_real,llm,name='u0')
    CALL allocate_field(f_q0,field_t,type_real,llm,nqtot,'q0')

    CALL allocate_field(f_dtheta_rhodz, field_t, type_real, llm, nqdyn, name='theta_rhodz0')
    CALL allocate_field(f_du,field_u,type_real,llm,name='u0')
    CALL allocate_field(f_dq,field_t,type_real,llm,nqtot,'q0')

    CALL allocate_field(f_rhodz, field_t, type_real, llm, name='rhodz')
    
    phys_smooth_tendency=.FALSE.
    CALL getin("phys_smooth_tendency",phys_smooth_tendency)
!$OMP END PARALLEL    
        
  END SUBROUTINE init_physics
  
  SUBROUTINE physics(it_,f_phis_, f_ps_, f_theta_rhodz_, f_u_, f_wflux_, f_q_)
  USE icosa
  USE field_mod
  USE mpipara
  USE omp_para
  !USE xios
  USE domain_mod
  USE time_mod
  USE disvert_mod
  IMPLICIT NONE
    INTEGER,INTENT(IN)    :: it_
    TYPE(t_field),POINTER :: f_phis_(:)
    TYPE(t_field),POINTER :: f_ps_(:)
    TYPE(t_field),POINTER :: f_theta_rhodz_(:)
    TYPE(t_field),POINTER :: f_u_(:)
    TYPE(t_field),POINTER :: f_wflux_(:)
    TYPE(t_field),POINTER :: f_q_(:)

    REAL(rstd),POINTER    :: theta_rhodz(:,:,:), theta_rhodz0(:,:,:), dtheta_rhodz(:,:,:) 
    REAL(rstd),POINTER    :: u(:,:), u0(:,:), du(:,:)
    REAL(rstd),POINTER    :: q(:,:,:),q0(:,:,:),dq(:,:,:)
    REAL(rstd),POINTER    :: ps(:)
    REAL(rstd),POINTER    :: rhodz(:,:)
    INTEGER :: ind, iq
    
    
!$OMP BARRIER
!$OMP MASTER
    f_phis=>f_phis_
    f_ps=>f_ps_
    f_theta_rhodz=>f_theta_rhodz_
    f_u=>f_u_
    f_wflux=>f_wflux_
    f_q=>f_q_
!$OMP END MASTER
!$OMP BARRIER

    IF (phys_smooth_tendency) THEN

      IF (MOD(it_,itau_physics)==1) THEN
        DO ind=1, ndomain
          IF (.NOT. assigned_domain(ind)) CYCLE
          CALL swap_dimensions(ind)
          CALL swap_geometry(ind)
          theta_rhodz=f_theta_rhodz(ind)
          theta_rhodz0=f_theta_rhodz0(ind)
          u=f_u(ind)
          u0=f_u0(ind)
          q=f_q(ind)
          q0=f_q0(ind)
          ps=f_ps(ind)
          rhodz=f_rhodz(ind)
       
          theta_rhodz0(:,ll_begin:ll_end,1)=theta_rhodz(:,ll_begin:ll_end,1)
          u0(:,ll_begin:ll_end)=u(:,ll_begin:ll_end)
          q0(:,ll_begin:ll_end,:)=q(:,ll_begin:ll_end,:)
          CALL compute_rhodz(.TRUE., ps, rhodz)
        ENDDO
     
!        IF (is_omp_master) CALL xios_timer_suspend("dynamico")
        it = it_-1 + itau_physics
        CALL external_physics
!        IF (is_omp_master) CALL xios_timer_resume("dynamico")

        DO ind=1, ndomain
          IF (.NOT. assigned_domain(ind)) CYCLE
          CALL swap_dimensions(ind) 
          CALL swap_geometry(ind) 
          theta_rhodz=f_theta_rhodz(ind)
          theta_rhodz0=f_theta_rhodz0(ind)
          u=f_u(ind)
          u0=f_u0(ind)
          q=f_q(ind)
          q0=f_q0(ind)
          dtheta_rhodz=f_dtheta_rhodz(ind)
          du=f_du(ind)
          dq=f_dq(ind)
          rhodz=f_rhodz(ind)
       
          dtheta_rhodz(:,ll_begin:ll_end,1)=(theta_rhodz(:,ll_begin:ll_end,1)-theta_rhodz0(:,ll_begin:ll_end,1))/itau_physics
          
          du(:,ll_begin:ll_end)=(u(:,ll_begin:ll_end)-u0(:,ll_begin:ll_end))/itau_physics
         
          DO iq=1, nqtot
            dq(:,ll_begin:ll_end,iq)=((q(:,ll_begin:ll_end,iq)-q0(:,ll_begin:ll_end,iq))/itau_physics)*rhodz(:,ll_begin:ll_end)
          ENDDO
         
          theta_rhodz(:,ll_begin:ll_end,1)=theta_rhodz0(:,ll_begin:ll_end,1)
          u(:,ll_begin:ll_end)=u0(:,ll_begin:ll_end)
          q(:,ll_begin:ll_end,:)=q0(:,ll_begin:ll_end,:)
        ENDDO
     ENDIF
   
     DO ind=1, ndomain
       IF (.NOT. assigned_domain(ind)) CYCLE
       CALL swap_dimensions(ind)
       CALL swap_geometry(ind)
         
       theta_rhodz=f_theta_rhodz(ind)
       u=f_u(ind)
       q=f_q(ind)
       dtheta_rhodz=f_dtheta_rhodz(ind)
       du=f_du(ind)
       dq=f_dq(ind)
       rhodz=f_rhodz(ind)
       ps=f_ps(ind)

       u(:,ll_begin:ll_end)=u(:,ll_begin:ll_end)+du(:,ll_begin:ll_end)
       theta_rhodz(:,ll_begin:ll_end,1)=theta_rhodz(:,ll_begin:ll_end,1)+dtheta_rhodz(:,ll_begin:ll_end,1)
       CALL compute_rhodz(.TRUE., ps, rhodz)
       DO iq=1, nqtot
         q(:,ll_begin:ll_end,iq)=q(:,ll_begin:ll_end,iq)+dq(:,ll_begin:ll_end,iq)/rhodz(:,ll_begin:ll_end)
       ENDDO
     ENDDO
!$OMP BARRIER
   
   ELSE
       
     IF (MOD(it_,itau_physics)==0) THEN
       it=it_
!       IF (is_omp_master) CALL xios_timer_suspend("dynamico")
       CALL external_physics
!       IF (is_omp_master) CALL xios_timer_resume("dynamico")
     ENDIF
   
   ENDIF      


  END SUBROUTINE physics
  
  
END MODULE physics_external_mod