source: LMDZ6/branches/Amaury_dev/libf/phylmd/ener_conserv.F90 @ 5185

SUBROUTINE ener_conserv(klon, klev, pdtphys, &
        puo, pvo, pto, qx, ivap, iliq, isol, &
        pun, pvn, ptn, pqn, pqln, pqsn, dtke, masse, exner, d_t_ec)

  !=============================================================
  ! Energy conservation
  ! Based on the TKE equation
  ! The M2 and N2 terms at the origin of TKE production are
  ! concerted into heating in the d_t_ec term
  ! Option 1 is the standard
  !        101 is for M2 term only
  !        101 for N2 term only
  !         -1 is a previours treatment for kinetic energy only
  !  FH (hourdin@lmd.jussieu.fr), 2013/04/25
  !=============================================================

  !=============================================================
  ! Declarations
  !=============================================================

  ! From module
  USE phys_local_var_mod, ONLY: d_u_vdf, d_v_vdf, d_t_vdf, d_u_ajs, d_v_ajs, d_t_ajs, &
          d_u_con, d_v_con, d_t_con, d_t_diss
  USE phys_local_var_mod, ONLY: d_t_eva, d_t_lsc, d_q_eva, d_q_lsc
  USE phys_local_var_mod, ONLY: d_u_oro, d_v_oro, d_u_lif, d_v_lif
  USE phys_local_var_mod, ONLY: du_gwd_hines, dv_gwd_hines, dv_gwd_front, dv_gwd_rando
  USE phys_state_var_mod, ONLY: du_gwd_front, du_gwd_rando
  USE phys_output_var_mod, ONLY: bils_ec, bils_ech, bils_tke, bils_kinetic, bils_enthalp, bils_latent, bils_diss
  USE add_phys_tend_mod, ONLY: fl_cor_ebil
  USE infotrac_phy, ONLY: nqtot
  USE lmdz_abort_physic, ONLY: abort_physic
  USE lmdz_clesphys
  USE lmdz_compbl, ONLY: iflag_pbl, iflag_pbl_split, iflag_order2_sollw, ifl_pbltree
  USE lmdz_yoethf
  USE lmdz_yomcst

  IMPLICIT NONE

  ! Arguments
  INTEGER, INTENT(IN) :: klon, klev
  REAL, INTENT(IN) :: pdtphys
  REAL, DIMENSION(klon, klev), INTENT(IN) :: puo, pvo, pto
  REAL, DIMENSION(klon, klev, nqtot), INTENT(IN) :: qx
  INTEGER, INTENT(IN) :: ivap, iliq, isol
  REAL, DIMENSION(klon, klev), INTENT(IN) :: pun, pvn, ptn, pqn, pqln, pqsn
  REAL, DIMENSION(klon, klev), INTENT(IN) :: masse, exner
  REAL, DIMENSION(klon, klev + 1), INTENT(IN) :: dtke

  REAL, DIMENSION(klon, klev), INTENT(OUT) :: d_t_ec

  ! Local
  INTEGER k, i
  REAL, DIMENSION(klon, klev + 1) :: fluxu, fluxv, fluxt
  REAL, DIMENSION(klon, klev + 1) :: dddu, dddv, dddt
  REAL, DIMENSION(klon, klev) :: d_u, d_v, d_t, zv, zu, d_t_ech, pqo, pql0, pqs0
  REAL ZRCPD

  CHARACTER*80 abort_message
  CHARACTER*20 :: modname

  modname = 'ener_conser'
  d_t_ec(:, :) = 0.

  IF(ivap == 0) CALL abort_physic (modname, 'can''t run without water vapour', 1)
  IF(iliq == 0) CALL abort_physic (modname, 'can''t run without liquid water', 1)
  pqo = qx(:, :, ivap)
  pql0 = qx(:, :, iliq)
  IF(isol /= 0) pqs0 = qx(:, :, isol)

  IF (iflag_ener_conserv==-1) THEN
    !+jld ec_conser
    DO k = 1, klev
      DO i = 1, klon
        IF (fl_cor_ebil > 0) THEN
          ZRCPD = RCPD * (1.0 + RVTMP2 * (pqn(i, k) + pqln(i, k) + pqsn(i, k)))
        ELSE
          ZRCPD = RCPD * (1.0 + RVTMP2 * pqn(i, k))
        ENDIF
        d_t_ec(i, k) = 0.5 / ZRCPD &
                * (puo(i, k)**2 + pvo(i, k)**2 - pun(i, k)**2 - pvn(i, k)**2)
      ENDDO
    ENDDO
    !-jld ec_conser

  ELSEIF (iflag_ener_conserv>=1) THEN

    IF (iflag_ener_conserv<=2) THEN
      !     PRINT*,'ener_conserv pbl=',iflag_pbl
      IF (iflag_pbl>=20 .AND. iflag_pbl<=27) THEN !d_t_diss accounts for conserv
        d_t(:, :) = d_t_ajs(:, :)   ! d_t_ajs = adjust + thermals
        d_u(:, :) = d_u_ajs(:, :) + d_u_con(:, :)
        d_v(:, :) = d_v_ajs(:, :) + d_v_con(:, :)
      ELSE
        d_t(:, :) = d_t_vdf(:, :) + d_t_ajs(:, :)   ! d_t_ajs = adjust + thermals
        d_u(:, :) = d_u_vdf(:, :) + d_u_ajs(:, :) + d_u_con(:, :)
        d_v(:, :) = d_v_vdf(:, :) + d_v_ajs(:, :) + d_v_con(:, :)
      ENDIF
    ELSEIF (iflag_ener_conserv==101) THEN
      d_t(:, :) = 0.
      d_u(:, :) = d_u_vdf(:, :) + d_u_ajs(:, :) + d_u_con(:, :)
      d_v(:, :) = d_v_vdf(:, :) + d_v_ajs(:, :) + d_v_con(:, :)
    ELSEIF (iflag_ener_conserv==110) THEN
      d_t(:, :) = d_t_vdf(:, :) + d_t_ajs(:, :)
      d_u(:, :) = 0.
      d_v(:, :) = 0.

    ELSEIF (iflag_ener_conserv==3) THEN
      d_t(:, :) = 0.
      d_u(:, :) = 0.
      d_v(:, :) = 0.
    ELSEIF (iflag_ener_conserv==4) THEN
      d_t(:, :) = 0.
      d_u(:, :) = d_u_vdf(:, :)
      d_v(:, :) = d_v_vdf(:, :)
    ELSEIF (iflag_ener_conserv==5) THEN
      d_t(:, :) = d_t_vdf(:, :)
      d_u(:, :) = d_u_vdf(:, :)
      d_v(:, :) = d_v_vdf(:, :)
    ELSEIF (iflag_ener_conserv==6) THEN
      d_t(:, :) = d_t_vdf(:, :)
      d_u(:, :) = d_u_vdf(:, :) + d_u_ajs(:, :)
      d_v(:, :) = d_v_vdf(:, :) + d_v_ajs(:, :)
    ELSEIF (iflag_ener_conserv==7) THEN
      d_t(:, :) = d_t_vdf(:, :) + d_t_ajs(:, :)
      d_u(:, :) = d_u_vdf(:, :) + d_u_ajs(:, :)
      d_v(:, :) = d_v_vdf(:, :) + d_v_ajs(:, :)
    ELSEIF (iflag_ener_conserv==8) THEN
      d_t(:, :) = d_t_vdf(:, :)
      d_u(:, :) = d_u_vdf(:, :) + d_u_ajs(:, :) + d_u_con(:, :)
      d_v(:, :) = d_v_vdf(:, :) + d_v_ajs(:, :) + d_v_con(:, :)
    ELSEIF (iflag_ener_conserv==9) THEN
      d_t(:, :) = d_t_vdf(:, :)
      d_u(:, :) = d_u_vdf(:, :) + d_u_ajs(:, :) + d_u_con(:, :) + d_u_oro(:, :)
      d_v(:, :) = d_v_vdf(:, :) + d_v_ajs(:, :) + d_v_con(:, :) + d_v_oro(:, :)
    ELSEIF (iflag_ener_conserv==10) THEN
      d_t(:, :) = d_t_vdf(:, :)
      d_u(:, :) = d_u_vdf(:, :) + d_u_ajs(:, :) + d_u_con(:, :) + d_u_oro(:, :) + d_u_lif(:, :)
      d_v(:, :) = d_v_vdf(:, :) + d_v_ajs(:, :) + d_v_con(:, :) + d_v_oro(:, :) + d_v_lif(:, :)
    ELSEIF (iflag_ener_conserv==11) THEN
      d_t(:, :) = d_t_vdf(:, :)
      d_u(:, :) = d_u_vdf(:, :) + d_u_ajs(:, :) + d_u_con(:, :) + d_u_oro(:, :) + d_u_lif(:, :)
      d_v(:, :) = d_v_vdf(:, :) + d_v_ajs(:, :) + d_v_con(:, :) + d_v_oro(:, :) + d_v_lif(:, :)
      IF (ok_hines) THEN
        d_u_vdf(:, :) = d_u_vdf(:, :) + du_gwd_hines(:, :)
        d_v_vdf(:, :) = d_v_vdf(:, :) + dv_gwd_hines(:, :)
      ENDIF
      IF (.NOT. ok_hines .AND. ok_gwd_rando) THEN
        d_u_vdf(:, :) = d_u_vdf(:, :) + du_gwd_front(:, :)
        d_v_vdf(:, :) = d_v_vdf(:, :) + dv_gwd_front(:, :)
      ENDIF
      IF (ok_gwd_rando) THEN
        d_u_vdf(:, :) = d_u_vdf(:, :) + du_gwd_rando(:, :)
        d_v_vdf(:, :) = d_v_vdf(:, :) + dv_gwd_rando(:, :)
      ENDIF
    ELSE
      abort_message = 'iflag_ener_conserv non prevu'
      CALL abort_physic (modname, abort_message, 1)
    ENDIF

    !----------------------------------------------------------------------------
    ! Two options wether we consider time integration in the energy conservation
    !----------------------------------------------------------------------------

    IF (iflag_ener_conserv==2) THEN
      zu(:, :) = puo(:, :)
      zv(:, :) = pvo(:, :)
    else
      IF (iflag_pbl>=20 .AND. iflag_pbl<=27) THEN
        zu(:, :) = puo(:, :) + d_u_vdf(:, :) + 0.5 * d_u(:, :)
        zv(:, :) = pvo(:, :) + d_v_vdf(:, :) + 0.5 * d_v(:, :)
      ELSE
        zu(:, :) = puo(:, :) + 0.5 * d_u(:, :)
        zv(:, :) = pvo(:, :) + 0.5 * d_v(:, :)
      ENDIF
    endif

    fluxu(:, klev + 1) = 0.
    fluxv(:, klev + 1) = 0.
    fluxt(:, klev + 1) = 0.

    DO k = klev, 1, -1
      fluxu(:, k) = fluxu(:, k + 1) + masse(:, k) * d_u(:, k)
      fluxv(:, k) = fluxv(:, k + 1) + masse(:, k) * d_v(:, k)
      fluxt(:, k) = fluxt(:, k + 1) + masse(:, k) * d_t(:, k) / exner(:, k)
    enddo

    dddu(:, 1) = 2 * zu(:, 1) * fluxu(:, 1)
    dddv(:, 1) = 2 * zv(:, 1) * fluxv(:, 1)
    dddt(:, 1) = (exner(:, 1) - 1.) * fluxt(:, 1)

    DO k = 2, klev
      dddu(:, k) = (zu(:, k) - zu(:, k - 1)) * fluxu(:, k)
      dddv(:, k) = (zv(:, k) - zv(:, k - 1)) * fluxv(:, k)
      dddt(:, k) = (exner(:, k) - exner(:, k - 1)) * fluxt(:, k)
    enddo
    dddu(:, klev + 1) = 0.
    dddv(:, klev + 1) = 0.
    dddt(:, klev + 1) = 0.

    DO k = 1, klev
      d_t_ech(:, k) = -(rcpd * (dddt(:, k) + dddt(:, k + 1))) / (2. * rcpd * masse(:, k))
      d_t_ec(:, k) = -(dddu(:, k) + dddu(:, k + 1) + dddv(:, k) + dddv(:, k + 1)) / (2. * rcpd * masse(:, k)) + d_t_ech(:, k)
    enddo

  ENDIF

  !================================================================
  !  Computation of integrated enthalpie and kinetic energy variation
  !  FH (hourdin@lmd.jussieu.fr), 2013/04/25
  !  bils_ec : energie conservation term
  !  bils_ech : part of this term linked to temperature
  !  bils_tke : change of TKE
  !  bils_diss : dissipation of TKE (when activated)
  !  bils_kinetic : change of kinetic energie of the column
  !  bils_enthalp : change of enthalpie
  !  bils_latent  : change of latent heat. Computed between
  !          after reevaporation (at the beginning of the physics)
  !          and before large scale condensation (fisrtilp)
  !================================================================

  bils_ec(:) = 0.
  bils_ech(:) = 0.
  bils_tke(:) = 0.
  bils_diss(:) = 0.
  bils_kinetic(:) = 0.
  bils_enthalp(:) = 0.
  bils_latent(:) = 0.
  DO k = 1, klev
    bils_ec(:) = bils_ec(:) - d_t_ec(:, k) * masse(:, k)
    bils_diss(:) = bils_diss(:) - d_t_diss(:, k) * masse(:, k)
    bils_kinetic(:) = bils_kinetic(:) + masse(:, k) * &
            (pun(:, k) * pun(:, k) + pvn(:, k) * pvn(:, k) &
                    - puo(:, k) * puo(:, k) - pvo(:, k) * pvo(:, k))
    bils_enthalp(:) = &
            bils_enthalp(:) + masse(:, k) * (ptn(:, k) - pto(:, k) + d_t_ec(:, k) - d_t_eva(:, k) - d_t_lsc(:, k))
    !    &  bils_enthalp(:)+masse(:,k)*(ptn(:,k)-pto(:,k)+d_t_ec(:,k))
    bils_latent(:) = bils_latent(:) + masse(:, k) * &
            !    &             (pqn(:,k)-pqo(:,k))
            (pqn(:, k) - pqo(:, k) - d_q_eva(:, k) - d_q_lsc(:, k))
  ENDDO
  bils_ec(:) = rcpd * bils_ec(:) / pdtphys
  bils_diss(:) = rcpd * bils_diss(:) / pdtphys
  bils_kinetic(:) = 0.5 * bils_kinetic(:) / pdtphys
  bils_enthalp(:) = rcpd * bils_enthalp(:) / pdtphys
  bils_latent(:) = rlvtt * bils_latent(:) / pdtphys
  !jyg<
  IF (iflag_pbl > 1) THEN
    DO k = 1, klev
      bils_tke(:) = bils_tke(:) + 0.5 * (dtke(:, k) + dtke(:, k + 1)) * masse(:, k)
    ENDDO
    bils_tke(:) = bils_tke(:) / pdtphys
  ENDIF  ! (iflag_pbl > 1)
  !>jyg

  IF (iflag_ener_conserv>=1) THEN
    bils_ech(:) = 0.
    DO k = 1, klev
      bils_ech(:) = bils_ech(:) - d_t_ech(:, k) * masse(:, k)
    ENDDO
    bils_ech(:) = rcpd * bils_ech(:) / pdtphys
  ENDIF

  RETURN

END