MODULE physics_mod

  USE field_mod

  PRIVATE

  INTEGER, PARAMETER :: phys_none=0, phys_HS94=1, phys_DCMIP=2

  INTEGER :: nb_extra_physics_2D, nb_extra_physics_3D, phys_type
  TYPE(t_field),POINTER :: f_extra_physics_2D(:), f_extra_physics_3D(:)

  CHARACTER(LEN=255) :: physics_type="automatic"
!$OMP THREADPRIVATE(physics_type)

  PUBLIC :: physics, init_physics
CONTAINS

  SUBROUTINE init_physics
    USE mpipara
    USE etat0_mod
    USE icosa
    USE physics_interface_mod
    USE physics_dcmip_mod, ONLY : init_physics_dcmip=>init_physics
    IMPLICIT NONE

    physics_type='automatic'
    CALL getin("physics",physics_type)

    SELECT CASE(TRIM(physics_type))
    CASE ('automatic')
       etat0_type='jablonowsky06'
       CALL getin("etat0",etat0_type)
       SELECT CASE(TRIM(etat0_type))
       CASE('held_suarez')
          phys_type = phys_HS94
       CASE DEFAULT
          IF(is_mpi_root) PRINT*,"NO PHYSICAL PACKAGE USED"
          phys_type = phys_none
       END SELECT

    CASE ('dcmip')
       CALL init_physics_dcmip
       nb_extra_physics_2D=1 ! precl
       nb_extra_physics_3D=0
       phys_type = phys_DCMIP
    CASE DEFAULT
       IF(is_mpi_root) PRINT*, 'init_physics : Bad selector for variable physics <',&
            TRIM(physics_type), '> options are <automatic>, <dcmip>, <dry>'
       STOP
    END SELECT

    IF(is_mpi_root) THEN
       PRINT *, 'phys_type = ',phys_type
       PRINT *, 'nb_extra_physics_2D = ', nb_extra_physics_2D
       PRINT *, 'nb_extra_physics_3D = ', nb_extra_physics_3D
    END IF
    IF(nb_extra_physics_2D>0) CALL allocate_field(f_extra_physics_2D,field_t,type_real,nb_extra_physics_2D)
    IF(nb_extra_physics_3D>0) CALL allocate_field(f_extra_physics_3D,field_t,type_real,llm,nb_extra_physics_3D)

  END SUBROUTINE init_physics

  SUBROUTINE physics(it,jD_cur,jH_cur,f_phis, f_ps, f_theta_rhodz, f_ue, f_q)
    USE icosa
    USE physics_interface_mod
 !   USE etat0_mod
    USE etat0_heldsz_mod
    IMPLICIT NONE
    INTEGER, INTENT(IN)   :: it
    REAL(rstd),INTENT(IN)::jD_cur,jH_cur
    TYPE(t_field),POINTER :: f_phis(:)
    TYPE(t_field),POINTER :: f_ps(:)
    TYPE(t_field),POINTER :: f_theta_rhodz(:)
    TYPE(t_field),POINTER :: f_ue(:)
    TYPE(t_field),POINTER :: f_q(:)
    REAL(rstd),POINTER :: phis(:)
    REAL(rstd),POINTER :: ps(:)
    REAL(rstd),POINTER :: theta_rhodz(:,:)
    REAL(rstd),POINTER :: ue(:,:)
    REAL(rstd),POINTER :: q(:,:,:)

    LOGICAL:: firstcall,lastcall
    INTEGER :: ind
    TYPE(physics_inout) :: args

    SELECT CASE(phys_type)
    CASE (phys_none)
       ! No physics, do nothing
    CASE(phys_HS94)
       CALL held_suarez(f_ps,f_theta_rhodz,f_ue) 
    CASE DEFAULT
       CALL transfert_request(f_ps,req_i1)
       CALL transfert_request(f_theta_rhodz,req_i1)
       CALL transfert_request(f_ue,req_e1_vect)
       CALL transfert_request(f_q,req_i1)
       
       args%dt_phys = dt
       DO ind=1,ndomain
          IF (.NOT. assigned_domain(ind)) CYCLE
          CALL swap_dimensions(ind)
          CALL swap_geometry(ind)
      
          phis=f_phis(ind)
          ps=f_ps(ind)
          theta_rhodz=f_theta_rhodz(ind)
          ue=f_ue(ind)
          q=f_q(ind)

          IF(nb_extra_physics_2D>0) args%extra_2D=f_extra_physics_2D(ind)
          IF(nb_extra_physics_3D>0) args%extra_3D=f_extra_physics_3D(ind)
          CALL physics_column(args, phis, ps, theta_rhodz, ue, q)
       ENDDO

       IF (mod(it,itau_out)==0 ) THEN
          IF(nb_extra_physics_2D>0) CALL writefield("extra_physics_2D",f_extra_physics_2D)
          IF(nb_extra_physics_3D>0) CALL writefield("extra_physics_3D",f_extra_physics_3D)
       ENDIF
    END SELECT

  END SUBROUTINE physics

  SUBROUTINE physics_column(args, phis, ps, theta_rhodz, ue, q)
    USE icosa
    USE wind_mod
    USE pression_mod
    USE theta2theta_rhodz_mod
    USE physics_interface_mod
    USE physics_dcmip_mod
    IMPLICIT NONE
    TYPE(physics_inout) :: args    
    REAL(rstd) :: phis(iim*jjm)
    REAL(rstd) :: ps(iim*jjm)
    REAL(rstd) :: theta_rhodz(iim*jjm,llm)
    REAL(rstd) :: ue(3*iim*jjm,llm)
    REAL(rstd), TARGET :: q(iim*jjm,llm,nqtot)
    ! local arrays
    REAL(rstd), TARGET :: lat(iim*jjm)
    REAL(rstd), TARGET :: lon(iim*jjm)
    REAL(rstd), TARGET :: p(iim*jjm,llm+1)
    REAL(rstd), TARGET :: Temp(iim*jjm,llm)
    REAL(rstd), TARGET :: ulon(iim*jjm,llm)
    REAL(rstd), TARGET :: ulat(iim*jjm,llm)
    REAL(rstd), TARGET :: dTemp(iim*jjm,llm)
    REAL(rstd), TARGET :: dulon(iim*jjm,llm)
    REAL(rstd), TARGET :: dulat(iim*jjm,llm)
    REAL(rstd), TARGET :: dq(iim*jjm,llm,nqtot)
    REAL(rstd) :: uc(iim*jjm,3,llm)  ! 3D velocity at cell centers

    INTEGER :: i,j,ij,l
    REAL(rstd) :: due, dt2

    DO j=jj_begin,jj_end
      DO i=ii_begin,ii_end
        ij=(j-1)*iim+i
        CALL xyz2lonlat(xyz_i(ij,:),lon(ij),lat(ij)) 
      ENDDO
    ENDDO

    ! Reconstruct wind vector at hexagons
    CALL compute_pression(ps,p,0)
    CALL compute_theta_rhodz2temperature(ps,theta_rhodz,Temp,0)
    CALL compute_wind_centered(ue,uc)
    CALL compute_wind_centered_lonlat_compound(uc, ulon, ulat)
    args%ngrid = iim*jjm
    args%lon => lon
    args%lat => lat
    args%p => p
    args%Temp => Temp
    args%ulon => ulon
    args%ulat => ulat
    args%q => q
    args%dTemp => dTemp
    args%dulon => dulon
    args%dulat => dulat
    args%dq => dq

    SELECT CASE(phys_type)
    CASE (phys_DCMIP)
       CALL compute_phys_wrap(args)
    END SELECT

    q = q + args%dt_phys * dq
    Temp = Temp + args%dt_phys * dTemp
    CALL compute_temperature2theta_rhodz(ps,Temp,theta_rhodz,0)
    
    ! Reconstruct wind tendencies at edges and add
    CALL compute_wind_centered_from_lonlat_compound(dulon,dulat,uc)
    dt2=.5*args%dt_phys
    DO l=1,llm
      DO j=jj_begin,jj_end
        DO i=ii_begin,ii_end
          ij=(j-1)*iim+i
          due = sum((uc(ij,:,l) + uc(ij+t_right,:,l))*ep_e(ij+u_right,:) )
          ue(ij+u_right,l) = ue(ij+u_right,l) + dt2*due

          due = sum( (uc(ij,:,l) + uc(ij+t_lup,:,l))*ep_e(ij+u_lup,:) )
          ue(ij+u_lup,l)=ue(ij+u_lup,l) + dt2*due

          due = sum( (uc(ij,:,l) + uc(ij+t_ldown,:,l))*ep_e(ij+u_ldown,:) )
          ue(ij+u_ldown,l)=ue(ij+u_ldown,l) + dt2*due
        ENDDO
      ENDDO
    ENDDO

  END SUBROUTINE physics_column

END MODULE physics_mod