! $Id: calwake.F90 1992 2014-03-05 13:19:12Z musat $

SUBROUTINE calwake(paprs, pplay, dtime, t, q, omgb, dt_dwn, dq_dwn, m_dwn, &
    m_up, dt_a, dq_a, sigd, wdt_pbl, wdq_pbl, udt_pbl, udq_pbl, wake_deltat, &
    wake_deltaq, wake_dth, wake_h, wake_s, wake_dens, wake_pe, wake_fip, &
    wake_gfl, dt_wake, dq_wake, wake_k, undi_t, undi_q, wake_omgbdth, &
    wake_dp_omgb, wake_dtke, wake_dqke, wake_dtpbl, wake_dqpbl, wake_omg, &
    wake_dp_deltomg, wake_spread, wake_cstar, wake_d_deltat_gw, wake_ddeltat, &
    wake_ddeltaq)
  ! **************************************************************
  ! *
  ! CALWAKE                                                     *
  ! interface avec le schema de calcul de la poche    *
  ! froide                                            *
  ! *
  ! written by   : CHERUY Frederique, 13/03/2000, 10.31.05      *
  ! modified by :  ROEHRIG Romain,    01/30/2007                *
  ! **************************************************************

  USE dimphy
  IMPLICIT NONE
  ! ======================================================================
  include "dimensions.h"
  ! #include "dimphy.h"
  include "YOMCST.h"

  ! Arguments
  ! ----------

  INTEGER i, l, ktopw(klon)
  REAL dtime

  REAL paprs(klon, klev+1), pplay(klon, klev)
  REAL t(klon, klev), q(klon, klev), omgb(klon, klev)
  REAL dt_dwn(klon, klev), dq_dwn(klon, klev), m_dwn(klon, klev)
  REAL m_up(klon, klev)
  REAL dt_a(klon, klev), dq_a(klon, klev)
  REAL wdt_pbl(klon, klev), wdq_pbl(klon, klev)
  REAL udt_pbl(klon, klev), udq_pbl(klon, klev)
  REAL wake_deltat(klon, klev), wake_deltaq(klon, klev)
  REAL dt_wake(klon, klev), dq_wake(klon, klev)
  REAL wake_d_deltat_gw(klon, klev)
  REAL wake_h(klon), wake_s(klon)
  REAL wake_dth(klon, klev)
  REAL wake_pe(klon), wake_fip(klon), wake_gfl(klon)
  REAL undi_t(klon, klev), undi_q(klon, klev)
  REAL wake_omgbdth(klon, klev), wake_dp_omgb(klon, klev)
  REAL wake_dtke(klon, klev), wake_dqke(klon, klev)
  REAL wake_dtpbl(klon, klev), wake_dqpbl(klon, klev)
  REAL wake_omg(klon, klev), wake_dp_deltomg(klon, klev)
  REAL wake_spread(klon, klev), wake_cstar(klon)
  REAL wake_ddeltat(klon, klev), wake_ddeltaq(klon, klev)
  REAL d_deltatw(klon, klev), d_deltaqw(klon, klev)
  INTEGER wake_k(klon)
  REAL sigd(klon)
  REAL wake_dens(klon)

  ! Variable internes
  ! -----------------

  REAL aire
  REAL p(klon, klev), ph(klon, klev+1), pi(klon, klev)
  REAL te(klon, klev), qe(klon, klev), omgbe(klon, klev+1)
  REAL dtdwn(klon, klev), dqdwn(klon, klev)
  REAL dta(klon, klev), dqa(klon, klev)
  REAL wdtpbl(klon, klev), wdqpbl(klon, klev)
  REAL udtpbl(klon, klev), udqpbl(klon, klev)
  REAL amdwn(klon, klev), amup(klon, klev)
  REAL dtw(klon, klev), dqw(klon, klev), dth(klon, klev)
  REAL d_deltat_gw(klon, klev)
  REAL dtls(klon, klev), dqls(klon, klev)
  REAL tu(klon, klev), qu(klon, klev)
  REAL hw(klon), sigmaw(klon), wape(klon), fip(klon), gfl(klon)
  REAL omgbdth(klon, klev+1), dp_omgb(klon, klev)
  REAL dtke(klon, klev), dqke(klon, klev)
  REAL dtpbl(klon, klev), dqpbl(klon, klev)
  REAL omg(klon, klev+1), dp_deltomg(klon, klev), spread(klon, klev)
  REAL cstar(klon)
  REAL sigd0(klon), wdens(klon)

  REAL rdcp

  ! print *, '-> calwake, wake_s ', wake_s(1)

  rdcp = 1./3.5


  ! -----------------------------------------------------------
  ! IM 290108     DO 999 i=1,klon   ! a vectoriser
  ! ----------------------------------------------------------


  DO l = 1, klev
    DO i = 1, klon
      p(i, l) = pplay(i, l)
      ph(i, l) = paprs(i, l)
      pi(i, l) = (pplay(i,l)/100000.)**rdcp

      te(i, l) = t(i, l)
      qe(i, l) = q(i, l)
      omgbe(i, l) = omgb(i, l)

      dtdwn(i, l) = dt_dwn(i, l)
      dqdwn(i, l) = dq_dwn(i, l)
      dta(i, l) = dt_a(i, l)
      dqa(i, l) = dq_a(i, l)
      wdtpbl(i, l) = wdt_pbl(i, l)
      wdqpbl(i, l) = wdq_pbl(i, l)
      udtpbl(i, l) = udt_pbl(i, l)
      udqpbl(i, l) = udq_pbl(i, l)
    END DO
  END DO

  omgbe(:, klev+1) = 0.

  DO i = 1, klon
    sigd0(i) = sigd(i)
  END DO
  ! print*, 'sigd0,sigd', sigd0, sigd(i)
  DO i = 1, klon
    ph(i, klev+1) = 0.
  END DO

  DO i = 1, klon
    ktopw(i) = wake_k(i)
  END DO

  DO l = 1, klev
    DO i = 1, klon
      dtw(i, l) = wake_deltat(i, l)
      dqw(i, l) = wake_deltaq(i, l)
    END DO
  END DO

  DO l = 1, klev
    DO i = 1, klon
      dtls(i, l) = dt_wake(i, l)
      dqls(i, l) = dq_wake(i, l)
    END DO
  END DO

  DO i = 1, klon
    hw(i) = wake_h(i)
    sigmaw(i) = wake_s(i)
  END DO

  ! fkc les flux de masses sont evalues aux niveaux et valent 0 a la surface
  ! fkc  on veut le flux de masse au milieu des couches

  DO l = 1, klev - 1
    DO i = 1, klon
      amdwn(i, l) = 0.5*(m_dwn(i,l)+m_dwn(i,l+1))
      amdwn(i, l) = (m_dwn(i,l+1))
    END DO
  END DO

  ! au sommet le flux de masse est nul

  DO i = 1, klon
    amdwn(i, klev) = 0.5*m_dwn(i, klev)
  END DO

  DO l = 1, klev
    DO i = 1, klon
      amup(i, l) = m_up(i, l)
    END DO
  END DO

  CALL wake(p, ph, pi, dtime, sigd0, te, qe, omgbe, dtdwn, dqdwn, amdwn, &
    amup, dta, dqa, wdtpbl, wdqpbl, udtpbl, udqpbl, dtw, dqw, dth, hw, &
    sigmaw, wape, fip, gfl, dtls, dqls, ktopw, omgbdth, dp_omgb, wdens, tu, &
    qu, dtke, dqke, dtpbl, dqpbl, omg, dp_deltomg, spread, cstar, &
    d_deltat_gw, d_deltatw, d_deltaqw)

  DO l = 1, klev
    DO i = 1, klon
      IF (ktopw(i)>0) THEN
        wake_deltat(i, l) = dtw(i, l)
        wake_deltaq(i, l) = dqw(i, l)
        wake_d_deltat_gw(i, l) = d_deltat_gw(i, l)
        wake_omgbdth(i, l) = omgbdth(i, l)
        wake_dp_omgb(i, l) = dp_omgb(i, l)
        wake_dtke(i, l) = dtke(i, l)
        wake_dqke(i, l) = dqke(i, l)
        wake_dtpbl(i, l) = dtpbl(i, l)
        wake_dqpbl(i, l) = dqpbl(i, l)
        wake_omg(i, l) = omg(i, l)
        wake_dp_deltomg(i, l) = dp_deltomg(i, l)
        wake_spread(i, l) = spread(i, l)
        wake_dth(i, l) = dth(i, l)
        dt_wake(i, l) = dtls(i, l)
        dq_wake(i, l) = dqls(i, l)
        undi_t(i, l) = tu(i, l)
        undi_q(i, l) = qu(i, l)
        wake_ddeltat(i, l) = d_deltatw(i, l)
        wake_ddeltaq(i, l) = d_deltaqw(i, l)
      ELSE
        wake_deltat(i, l) = 0.
        wake_deltaq(i, l) = 0.
        wake_d_deltat_gw(i, l) = 0.
        wake_omgbdth(i, l) = 0.
        wake_dp_omgb(i, l) = 0.
        wake_dtke(i, l) = 0.
        wake_dqke(i, l) = 0.
        wake_dtpbl(i, l) = 0.
        wake_dqpbl(i, l) = 0.
        wake_omg(i, l) = 0.
        wake_dp_deltomg(i, l) = 0.
        wake_spread(i, l) = 0.
        wake_dth(i, l) = 0.
        dt_wake(i, l) = 0.
        dq_wake(i, l) = 0.
        undi_t(i, l) = te(i, l)
        undi_q(i, l) = qe(i, l)
        wake_ddeltat(i, l) = 0.
        wake_ddeltaq(i, l) = 0.
      END IF
    END DO
  END DO

  DO i = 1, klon
    wake_h(i) = hw(i)
    wake_s(i) = sigmaw(i)
    wake_pe(i) = wape(i)
    wake_fip(i) = fip(i)
    wake_gfl(i) = gfl(i)
    wake_k(i) = ktopw(i)
    wake_cstar(i) = cstar(i)
    wake_dens(i) = wdens(i)
  END DO

  RETURN
END SUBROUTINE calwake

SUBROUTINE calwake_scal(paprs, pplay, dtime, t, q, omgb, dt_dwn, dq_dwn, &
    m_dwn, m_up, dt_a, dq_a, sigd, wdt_pbl, wdq_pbl, udt_pbl, udq_pbl, &
    wake_deltat, wake_deltaq, wake_dth, wake_h, wake_s, wake_dens, wake_pe, &
    wake_fip, wake_gfl, dt_wake, dq_wake, wake_k, undi_t, undi_q, &
    wake_omgbdth, wake_dp_omgb, wake_dtke, wake_dqke, wake_dtpbl, wake_dqpbl, &
    wake_omg, wake_dp_deltomg, wake_spread, wake_cstar, wake_d_deltat_gw, &
    wake_ddeltat, wake_ddeltaq)
  ! **************************************************************
  ! *
  ! CALWAKE                                                     *
  ! interface avec le schema de calcul de la poche    *
  ! froide                                            *
  ! *
  ! written by   : CHERUY Frederique, 13/03/2000, 10.31.05      *
  ! modified by :  ROEHRIG Romain,    01/30/2007                *
  ! **************************************************************

  USE dimphy
  IMPLICIT NONE
  ! ======================================================================

  include "dimensions.h"
  ! ccc#include "dimphy.h"
  include "YOMCST.h"

  ! Arguments
  ! ----------

  INTEGER i, l, ktopw
  REAL dtime

  REAL paprs(klon, klev+1), pplay(klon, klev)
  REAL t(klon, klev), q(klon, klev), omgb(klon, klev)
  REAL dt_dwn(klon, klev), dq_dwn(klon, klev), m_dwn(klon, klev)
  REAL m_up(klon, klev)
  REAL dt_a(klon, klev), dq_a(klon, klev)
  REAL wdt_pbl(klon, klev), wdq_pbl(klon, klev)
  REAL udt_pbl(klon, klev), udq_pbl(klon, klev)
  REAL wake_deltat(klon, klev), wake_deltaq(klon, klev)
  REAL dt_wake(klon, klev), dq_wake(klon, klev)
  REAL wake_d_deltat_gw(klon, klev)
  REAL wake_h(klon), wake_s(klon)
  REAL wake_dth(klon, klev)
  REAL wake_pe(klon), wake_fip(klon), wake_gfl(klon)
  REAL undi_t(klon, klev), undi_q(klon, klev)
  REAL wake_omgbdth(klon, klev), wake_dp_omgb(klon, klev)
  REAL wake_dtke(klon, klev), wake_dqke(klon, klev)
  REAL wake_dtpbl(klon, klev), wake_dqpbl(klon, klev)
  REAL wake_omg(klon, klev+1), wake_dp_deltomg(klon, klev)
  REAL wake_spread(klon, klev), wake_cstar(klon)
  REAL wake_ddeltat(klon, klev), wake_ddeltaq(klon, klev)
  REAL d_deltatw(klev), d_deltaqw(klev)
  INTEGER wake_k(klon)
  REAL sigd(klon)
  REAL wake_dens(klon)

  ! Variable internes
  ! -----------------

  REAL aire
  REAL p(klev), ph(klev+1), pi(klev)
  REAL te(klev), qe(klev), omgbe(klev), dtdwn(klev), dqdwn(klev)
  REAL dta(klev), dqa(klev)
  REAL wdtpbl(klev), wdqpbl(klev)
  REAL udtpbl(klev), udqpbl(klev)
  REAL amdwn(klev), amup(klev)
  REAL dtw(klev), dqw(klev), dth(klev), d_deltat_gw(klev)
  REAL dtls(klev), dqls(klev)
  REAL tu(klev), qu(klev)
  REAL hw, sigmaw, wape, fip, gfl
  REAL omgbdth(klev), dp_omgb(klev)
  REAL dtke(klev), dqke(klev)
  REAL dtpbl(klev), dqpbl(klev)
  REAL omg(klev+1), dp_deltomg(klev), spread(klev), cstar
  REAL sigd0, wdens

  REAL rdcp

  ! print *, '-> calwake, wake_s ', wake_s(1)

  rdcp = 1./3.5

  ! -----------------------------------------------------------
  DO i = 1, klon ! a vectoriser
    ! ----------------------------------------------------------


    DO l = 1, klev
      p(l) = pplay(i, l)
      ph(l) = paprs(i, l)
      pi(l) = (pplay(i,l)/100000.)**rdcp

      te(l) = t(i, l)
      qe(l) = q(i, l)
      omgbe(l) = omgb(i, l)

      dtdwn(l) = dt_dwn(i, l)
      dqdwn(l) = dq_dwn(i, l)
      dta(l) = dt_a(i, l)
      dqa(l) = dq_a(i, l)
      wdtpbl(l) = wdt_pbl(i, l)
      wdqpbl(l) = wdq_pbl(i, l)
      udtpbl(l) = udt_pbl(i, l)
      udqpbl(l) = udq_pbl(i, l)
    END DO

    sigd0 = sigd(i)
    ! print*, 'sigd0,sigd', sigd0, sigd(i)
    ph(klev+1) = 0.

    ktopw = wake_k(i)

    DO l = 1, klev
      dtw(l) = wake_deltat(i, l)
      dqw(l) = wake_deltaq(i, l)
    END DO

    DO l = 1, klev
      dtls(l) = dt_wake(i, l)
      dqls(l) = dq_wake(i, l)
    END DO

    hw = wake_h(i)
    sigmaw = wake_s(i)

    ! fkc les flux de masses sont evalues aux niveaux et valent 0 a la
    ! surface
    ! fkc  on veut le flux de masse au milieu des couches

    DO l = 1, klev - 1
      amdwn(l) = 0.5*(m_dwn(i,l)+m_dwn(i,l+1))
      amdwn(l) = (m_dwn(i,l+1))
    END DO

    ! au sommet le flux de masse est nul

    amdwn(klev) = 0.5*m_dwn(i, klev)

    DO l = 1, klev
      amup(l) = m_up(i, l)
    END DO

    CALL wake_scal(p, ph, pi, dtime, sigd0, te, qe, omgbe, dtdwn, dqdwn, &
      amdwn, amup, dta, dqa, wdtpbl, wdqpbl, udtpbl, udqpbl, dtw, dqw, dth, &
      hw, sigmaw, wape, fip, gfl, dtls, dqls, ktopw, omgbdth, dp_omgb, wdens, &
      tu, qu, dtke, dqke, dtpbl, dqpbl, omg, dp_deltomg, spread, cstar, &
      d_deltat_gw, d_deltatw, d_deltaqw)

    IF (ktopw>0) THEN
      DO l = 1, klev
        wake_deltat(i, l) = dtw(l)
        wake_deltaq(i, l) = dqw(l)
        wake_d_deltat_gw(i, l) = d_deltat_gw(l)
        wake_omgbdth(i, l) = omgbdth(l)
        wake_dp_omgb(i, l) = dp_omgb(l)
        wake_dtke(i, l) = dtke(l)
        wake_dqke(i, l) = dqke(l)
        wake_dtpbl(i, l) = dtpbl(l)
        wake_dqpbl(i, l) = dqpbl(l)
        wake_omg(i, l) = omg(l)
        wake_dp_deltomg(i, l) = dp_deltomg(l)
        wake_spread(i, l) = spread(l)
        wake_dth(i, l) = dth(l)
        dt_wake(i, l) = dtls(l)
        dq_wake(i, l) = dqls(l)
        undi_t(i, l) = tu(l)
        undi_q(i, l) = qu(l)
        wake_ddeltat(i, l) = d_deltatw(l)
        wake_ddeltaq(i, l) = d_deltaqw(l)
      END DO
    ELSE
      DO l = 1, klev
        wake_deltat(i, l) = 0.
        wake_deltaq(i, l) = 0.
        wake_d_deltat_gw(i, l) = 0.
        wake_omgbdth(i, l) = 0.
        wake_dp_omgb(i, l) = 0.
        wake_dtke(i, l) = 0.
        wake_dqke(i, l) = 0.
        wake_omg(i, l) = 0.
        wake_dp_deltomg(i, l) = 0.
        wake_spread(i, l) = 0.
        wake_dth(i, l) = 0.
        dt_wake(i, l) = 0.
        dq_wake(i, l) = 0.
        undi_t(i, l) = te(l)
        undi_q(i, l) = qe(l)
      END DO
    END IF

    wake_h(i) = hw
    wake_s(i) = sigmaw
    wake_pe(i) = wape
    wake_fip(i) = fip
    wake_gfl(i) = gfl
    wake_k(i) = ktopw
    wake_cstar(i) = cstar
    wake_dens(i) = wdens

  END DO

  RETURN
END SUBROUTINE calwake_scal