MODULE physics_dcmip2016_mod
   USE ICOSA
   PRIVATE

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

   TYPE(t_field), SAVE, POINTER  :: f_precl(:)
   REAL(rstd), SAVE, ALLOCATABLE :: precl_packed(:)
!$OMP THREADPRIVATE(precl_packed)

   TYPE(t_field), SAVE, POINTER  :: f_Q1(:)
   TYPE(t_field), SAVE, POINTER  :: f_Q2(:)
   TYPE(t_field), SAVE, POINTER  :: f_PS(:)
   TYPE(t_field), SAVE, POINTER  :: f_rhodz(:)
   TYPE(t_field), SAVE, POINTER  :: f_Q1_col_int(:)
   TYPE(t_field), SAVE, POINTER  :: f_Q2_col_int(:)
   PUBLIC :: init_physics, full_physics, write_physics

   INTEGER, PARAMETER :: dry_baroclinic = 0
   INTEGER, PARAMETER :: moist_baroclinic_full = 1
   INTEGER, PARAMETER :: moist_baroclinic_kessler = 2
   INTEGER, PARAMETER :: cyclone = 3
   INTEGER, PARAMETER :: supercell = 4

   LOGICAL, SAVE :: PBL   !  boundary layer
   !  True : George Bryan
   !  False : Reed & Jablonowsi
   !$OMP THREADPRIVATE(PBL)
CONTAINS

   SUBROUTINE init_physics
      USE physics_interface_mod
      USE omp_para
      IMPLICIT NONE
      INTEGER :: ngrid
      CHARACTER(LEN=255) :: testcase_str

      testcase_str = 'undefined'
      CALL getin("physics_dcmip2016", testcase_str)

      SELECT CASE (TRIM(testcase_str))
      CASE ('dry_baroclinic')
         testcase = dry_baroclinic
      CASE ('moist_baroclinic_full')
         testcase = moist_baroclinic_full
      CASE ('moist_baroclinic_kessler')
         testcase = moist_baroclinic_kessler
      CASE ('cyclone')
         testcase = cyclone
      CASE ('supercell')
         testcase = supercell
      CASE DEFAULT
         PRINT *, 'Bad selector for dcmip2016 test case <', testcase_str, &
            '> options are <dry_baroclinic>, <moist_baroclinic_full>, <moist_baroclinic_kessler>, <cyclone>, <supercell>'
         STOP
      END SELECT

      PBL = .FALSE.
      CALL getin("physics_dcmip2016_PBL", PBL)
      IF (is_master) THEN
         IF (PBL) THEN
            PRINT *, 'PBL=.TRUE., Bryan PBL activated'
         ELSE
            PRINT *, 'PBL=.FALSE., Reed & Jablonowski PBL activated'
         END IF
      END IF

      ! Physics-specific data
      ngrid = physics_inout%ngrid
      ALLOCATE (precl_packed(ngrid))
      CALL allocate_field(f_precl, field_t, type_real)
      CALL allocate_field(f_Q1, field_t, type_real, llm)
      CALL allocate_field(f_Q2, field_t, type_real, llm)
      CALL allocate_field(f_PS, field_t, type_real)
      CALL allocate_field(f_rhodz, field_t, type_real, llm)
      CALL allocate_field(f_Q1_col_int, field_t, type_real)
      CALL allocate_field(f_Q2_col_int, field_t, type_real)

      PRINT *, 'init_physics_new', SIZE(physics_inout%Ai)
   END SUBROUTINE init_physics

   SUBROUTINE full_physics
      USE physics_interface_mod
      CALL compute_physics(physics_inout%ngrid, physics_inout%dt_phys, &
                           physics_inout%lon, physics_inout%lat, physics_inout%p, physics_inout%pk, physics_inout%Temp, &
                           physics_inout%ulon, physics_inout%ulat, physics_inout%q(:, :, 1:5), &
                           physics_inout%dTemp, physics_inout%dulon, physics_inout%dulat, &
                           physics_inout%dq(:, :, 1:5), precl_packed)
   END SUBROUTINE full_physics

   SUBROUTINE write_physics
      USE output_field_mod
      USE physics_interface_mod
      use disvert_mod
      REAL(rstd), POINTER :: Q1(:, :)
      REAL(rstd), POINTER :: Q2(:, :)
      REAL(rstd), POINTER :: PS(:)
      REAL(rstd), POINTER :: rhodz(:, :)
      REAL(rstd), POINTER :: Q1_col_int(:)
      REAL(rstd), POINTER :: Q2_col_int(:)

      CALL unpack_field(f_precl, precl_packed)
      CALL output_field("precl", f_precl)

      CALL unpack_field(f_Q1, physics_inout%q(:, :, 4))
      CALL unpack_field(f_Q2, physics_inout%q(:, :, 5))
      CALL unpack_field(f_ps, physics_inout%p(:, 1))

      DO ind = 1, ndomain
         IF (.NOT. assigned_domain(ind)) CYCLE
         Q1 = f_Q1(ind)
         Q2 = f_Q2(ind)
         Q1_col_int = f_Q1_col_int(ind)
         Q2_col_int = f_Q2_col_int(ind)
         PS = f_PS(ind)
         rhodz = f_rhodz(ind)
         DO l = 1, llm
            rhodz(:, l) = (ap(l) - ap(l + 1) + (bp(l) - bp(l + 1))*ps(:))/g
         END DO
         Q1_col_int = SUM(Q1*rhodz, 2)/SUM(rhodz, 2)
         Q2_col_int = SUM(Q2*rhodz, 2)/SUM(rhodz, 2)
      END DO

      CALL output_field("Q1_col_int", f_Q1_col_int)
      CALL output_field("Q2_col_int", f_Q2_col_int)

   END SUBROUTINE write_physics

   SUBROUTINE compute_physics(ngrid, dt_phys, lon, lat, p, pk, Temp, u, v, q, dTemp, du, dv, dq, precl)
      USE icosa
      USE dcmip2016_simple_physics_mod
      USE dcmip2016_kessler_physic_mod
      USE earth_const
      USE terminator
      IMPLICIT NONE
      INTEGER    :: ngrid
      REAL(rstd) :: lat(ngrid)
      REAL(rstd) :: lon(ngrid)
      REAL(rstd) :: ps(ngrid)
      REAL(rstd) :: precl(ngrid)
      ! arguments with bottom-up indexing (DYNAMICO)
      REAL(rstd) :: p(ngrid, llm + 1)
      REAL(rstd) :: pk(ngrid, llm)
      REAL(rstd) :: Temp(ngrid, llm)
      REAL(rstd) :: u(ngrid, llm)
      REAL(rstd) :: v(ngrid, llm)
      REAL(rstd) :: q(ngrid, llm, 5)
      REAL(rstd) :: dTemp(ngrid, llm)
      REAL(rstd) :: du(ngrid, llm)
      REAL(rstd) :: dv(ngrid, llm)
      REAL(rstd) :: dq(ngrid, llm, 5)
      ! local arrays with top-down vertical indexing (DCMIP)
      REAL(rstd) :: pint(ngrid, llm + 1)
      REAL(rstd) :: pmid(ngrid, llm)
      REAL(rstd) :: pdel(ngrid, llm)
      REAL(rstd) :: Tfi(ngrid, llm)
      REAL(rstd) :: ufi(ngrid, llm)
      REAL(rstd) :: vfi(ngrid, llm)
      REAL(rstd) :: qfi(ngrid, llm, 5)

      REAL(rstd)  :: rho(llm), z(llm), theta(llm), qv(llm), qc(llm), qr(llm)
      REAL(rstd)  :: lastz
      REAL(rstd)  :: dcl1, dcl2
      INTEGER :: l, ll, ij
      REAL(rstd) :: dt_phys, inv_dt
      INTEGER :: simple_physic_testcase

      ! prepare input fields and mirror vertical index
      ps(:) = p(:, 1) ! surface pressure

      DO l = 1, llm + 1
         DO ij = 1, ngrid
            pint(ij, l) = p(ij, llm + 2 - l)
         END DO
      END DO

      DO l = 1, llm
         ll = llm + 1 - l
         DO ij = 1, ngrid
            pmid(ij, l) = 0.5*(pint(ij, l) + pint(ij, l + 1)) ! Pressure inside layers
            pdel(ij, l) = pint(ij, l + 1) - pint(ij, l)       ! Pressure difference between two layers
            ufi(ij, l) = u(ij, ll)
            vfi(ij, l) = v(ij, ll)
            qfi(ij, l, :) = q(ij, ll, :)
            IF (physics_thermo == thermo_fake_moist) THEN
               Tfi(ij, l) = Temp(ij, ll)/(1 + 0.608*qfi(ij, l, 1))
            ELSE
               Tfi(ij, l) = Temp(ij, ll)
            END IF
         END DO
      END DO

      precl = 0.
      IF (testcase == moist_baroclinic_full .OR. testcase == cyclone) THEN
         IF (testcase == moist_baroclinic_full) simple_physic_testcase = 1
         IF (testcase == cyclone) simple_physic_testcase = 0
         CALL simple_physics(ngrid, llm, dt_phys, lat, tfi, qfi(:, :, 1), ufi, vfi, pmid, pint, pdel, 1./pdel, ps, precl, &
                             simple_physic_testcase, .FALSE., PBL)
      END IF

      IF (testcase == moist_baroclinic_full .OR. testcase == moist_baroclinic_kessler &
      &   .OR. testcase == cyclone .OR. testcase == supercell) THEN
         DO ij = 1, ngrid
            lastz = 0
            DO l = 1, llm
               ll = llm + 1 - l
               rho(l) = pmid(ij, ll)/(287*Tfi(ij, ll)*(1 + 0.608*qfi(ij, ll, 1)))
               z(l) = lastz
               lastz = lastz + (p(ij, l) - p(ij, l + 1))/g/rho(l)
               theta(l) = Tfi(ij, ll)/(pk(ij, l)/cpp)
            END DO

            DO l = 1, llm - 1
               z(l) = 0.5*(z(l) + z(l + 1))
            END DO
            z(llm) = z(llm) + (z(llm) - z(llm - 1))

            qv(:) = max(qfi(ij, llm:1:-1, 1), 0.)
            qc(:) = max(qfi(ij, llm:1:-1, 2), 0.)
            qr(:) = max(qfi(ij, llm:1:-1, 3), 0.)

            CALL KESSLER(theta(:), qv, qc, qr, rho(:), &
                         pk(ij, :)/cpp, dt_phys, z(:), llm, precl(ij))

            DO l = 1, llm
               ll = llm + 1 - l
               Tfi(ij, ll) = theta(l)*(pk(ij, l)/cpp)
            END DO

            qfi(ij, :, 1) = qv(llm:1:-1)
            qfi(ij, :, 2) = qc(llm:1:-1)
            qfi(ij, :, 3) = qr(llm:1:-1)

         END DO
      END IF

      DO l = 1, llm
         ll = llm + 1 - l
         DO ij = 1, ngrid
            CALL tendency_terminator(lat(ij), lon(ij), qfi(ij, ll, 4), qfi(ij, ll, 5), dt_phys, dcl1, dcl2)
            qfi(ij, ll, 4) = qfi(ij, ll, 4) + dt_phys*dcl1
            qfi(ij, ll, 5) = qfi(ij, ll, 5) + dt_phys*dcl2
         END DO
      END DO

      ! Mirror vertical index and compute tendencies
      inv_dt = 1./dt_phys
      DO l = 1, llm
         ll = llm + 1 - l
         DO ij = 1, ngrid
            IF (physics_thermo == thermo_fake_moist) THEN
               dTemp(ij, l) = inv_dt*(Tfi(ij, ll)*(1 + 0.608*qfi(ij, ll, 1)) - Temp(ij, l))
            ELSE
               dTemp(ij, l) = inv_dt*(Tfi(ij, ll) - Temp(ij, l))
            END IF

            du(ij, l) = inv_dt*(ufi(ij, ll) - u(ij, l))
            dv(ij, l) = inv_dt*(vfi(ij, ll) - v(ij, l))
            dq(ij, l, :) = inv_dt*(qfi(ij, ll, :) - q(ij, l, :))
         END DO
      END DO

   END SUBROUTINE compute_physics

END MODULE physics_dcmip2016_mod
