Head

GCC Code Coverage Report

Directory:	./		Exec	Total	Coverage
File:	phylmd/climb_hq_mod.F90	Lines:	101	147	68.7 %
Date:	2023-06-30 12:51:15	Branches:	187	358	52.2 %


MODULE climb_hq_mod
!
! Module to solve the verctical diffusion of "q" and "H";
! specific humidity and potential energi.
!
  USE dimphy

  IMPLICIT NONE
  SAVE
  PRIVATE
  PUBLIC :: climb_hq_down, climb_hq_up, d_h_col_vdf, f_h_bnd

  REAL, DIMENSION(:,:), ALLOCATABLE :: gamaq, gamah
  !$OMP THREADPRIVATE(gamaq,gamah)
  REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_Q, Dcoef_Q
  !$OMP THREADPRIVATE(Ccoef_Q, Dcoef_Q)
  REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_H, Dcoef_H
  !$OMP THREADPRIVATE(Ccoef_H, Dcoef_H)
  REAL, DIMENSION(:), ALLOCATABLE   :: Acoef_Q, Bcoef_Q
  !$OMP THREADPRIVATE(Acoef_Q, Bcoef_Q)
  REAL, DIMENSION(:), ALLOCATABLE   :: Acoef_H, Bcoef_H
  !$OMP THREADPRIVATE(Acoef_H, Bcoef_H)
  REAL, DIMENSION(:,:), ALLOCATABLE :: Kcoefhq
  !$OMP THREADPRIVATE(Kcoefhq)
  REAL, SAVE, DIMENSION(:,:), ALLOCATABLE :: h_old ! for diagnostics, h before solving diffusion
  !$OMP THREADPRIVATE(h_old)
  REAL, SAVE, DIMENSION(:), ALLOCATABLE :: d_h_col_vdf ! for diagnostics, vertical integral of enthalpy change
  !$OMP THREADPRIVATE(d_h_col_vdf)
  REAL, SAVE, DIMENSION(:), ALLOCATABLE :: f_h_bnd ! for diagnostics, enthalpy flux at surface
  !$OMP THREADPRIVATE(f_h_bnd)

CONTAINS
!
!****************************************************************************************
!
  SUBROUTINE climb_hq_down(knon, coefhq, paprs, pplay, &
       delp, temp, q, dtime, &
!!! nrlmd le 02/05/2011
       Ccoef_H_out, Ccoef_Q_out, Dcoef_H_out, Dcoef_Q_out, &
       Kcoef_hq_out, gama_q_out, gama_h_out, &
!!!
       Acoef_H_out, Acoef_Q_out, Bcoef_H_out, Bcoef_Q_out)

! This routine calculates recursivly the coefficients C and D
! for the quantity X=[Q,H] in equation X(k) = C(k) + D(k)*X(k-1), where k is
! the index of the vertical layer.
!
! Input arguments
!****************************************************************************************
    INTEGER, INTENT(IN)                      :: knon
    REAL, DIMENSION(klon,klev), INTENT(IN)   :: coefhq
    REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay
    REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
    REAL, DIMENSION(klon,klev), INTENT(IN)   :: temp, delp  ! temperature
    REAL, DIMENSION(klon,klev), INTENT(IN)   :: q
    REAL, INTENT(IN)                         :: dtime

! Output arguments
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef_H_out
    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef_Q_out
    REAL, DIMENSION(klon), INTENT(OUT)       :: Bcoef_H_out
    REAL, DIMENSION(klon), INTENT(OUT)       :: Bcoef_Q_out

!!! nrlmd le 02/05/2011
    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef_H_out
    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef_Q_out
    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Dcoef_H_out
    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Dcoef_Q_out
    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Kcoef_hq_out
    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: gama_q_out
    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: gama_h_out
!!!

! Local variables
!****************************************************************************************
    LOGICAL, SAVE                            :: first=.TRUE.
    !$OMP THREADPRIVATE(first)
! JLD now renamed h_old and declared in module
!    REAL, DIMENSION(klon,klev)               :: local_H
    REAL, DIMENSION(klon)                    :: psref
    REAL                                     :: delz, pkh
    INTEGER                                  :: k, i, ierr
! Include
!****************************************************************************************
    INCLUDE "YOMCST.h"
    INCLUDE "compbl.h"


!****************************************************************************************
! 1)
! Allocation at first time step only
!
!****************************************************************************************

    IF (first) THEN
       first=.FALSE.
       ALLOCATE(Ccoef_Q(klon,klev), STAT=ierr)
       IF ( ierr /= 0 )  PRINT*,' pb in allloc Ccoef_Q, ierr=', ierr

       ALLOCATE(Dcoef_Q(klon,klev), STAT=ierr)
       IF ( ierr /= 0 )  PRINT*,' pb in allloc Dcoef_Q, ierr=', ierr

       ALLOCATE(Ccoef_H(klon,klev), STAT=ierr)
       IF ( ierr /= 0 )  PRINT*,' pb in allloc Ccoef_H, ierr=', ierr

       ALLOCATE(Dcoef_H(klon,klev), STAT=ierr)
       IF ( ierr /= 0 )  PRINT*,' pb in allloc Dcoef_H, ierr=', ierr

       ALLOCATE(Acoef_Q(klon), Bcoef_Q(klon), Acoef_H(klon), Bcoef_H(klon), STAT=ierr)
       IF ( ierr /= 0 )  PRINT*,' pb in allloc Acoef_X and Bcoef_X, ierr=', ierr

       ALLOCATE(Kcoefhq(klon,klev), STAT=ierr)
       IF ( ierr /= 0 )  PRINT*,' pb in allloc Kcoefhq, ierr=', ierr

       ALLOCATE(gamaq(1:klon,2:klev), STAT=ierr)
       IF ( ierr /= 0 ) PRINT*,' pb in allloc gamaq, ierr=', ierr

       ALLOCATE(gamah(1:klon,2:klev), STAT=ierr)
       IF ( ierr /= 0 ) PRINT*,' pb in allloc gamah, ierr=', ierr

       ALLOCATE(h_old(klon,klev), STAT=ierr)
       IF ( ierr /= 0 )  PRINT*,' pb in allloc h_old, ierr=', ierr

       ALLOCATE(d_h_col_vdf(klon), STAT=ierr)
       IF ( ierr /= 0 )  PRINT*,' pb in allloc d_h_col_vdf, ierr=', ierr

       ALLOCATE(f_h_bnd(klon), STAT=ierr)
       IF ( ierr /= 0 )  PRINT*,' pb in allloc f_h_bnd, ierr=', ierr
    END IF

!****************************************************************************************
! 2)
! Definition of the coeficient K
!
!****************************************************************************************
    Kcoefhq(:,:) = 0.0
    DO k = 2, klev
       DO i = 1, knon
          Kcoefhq(i,k) = &
               coefhq(i,k)*RG*RG*dtime /(pplay(i,k-1)-pplay(i,k)) &
               *(paprs(i,k)*2/(temp(i,k)+temp(i,k-1))/RD)**2
       ENDDO
    ENDDO

!****************************************************************************************
! 3)
! Calculation of gama for "Q" and "H"
!
!****************************************************************************************
!   surface pressure is used as reference
    psref(:) = paprs(:,1)

!   definition of gama
    IF (iflag_pbl == 1) THEN
       gamaq(:,:) = 0.0
       gamah(:,:) = -1.0e-03
       gamah(:,2) = -2.5e-03

! conversion de gama
       DO k = 2, klev
          DO i = 1, knon
             delz = RD * (temp(i,k-1)+temp(i,k)) / &
                    2.0 / RG / paprs(i,k) * (pplay(i,k-1)-pplay(i,k))
             pkh  = (psref(i)/paprs(i,k))**RKAPPA

! convertie gradient verticale d'humidite specifique en difference d'humidite specifique entre centre de couches
             gamaq(i,k) = gamaq(i,k) * delz
! convertie gradient verticale de temperature en difference de temperature potentielle entre centre de couches
             gamah(i,k) = gamah(i,k) * delz * RCPD * pkh
          ENDDO
       ENDDO

    ELSE
       gamaq(:,:) = 0.0
       gamah(:,:) = 0.0
    ENDIF


!****************************************************************************************
! 4)
! Calculte the coefficients C and D for specific humidity, q
!
!****************************************************************************************

    CALL calc_coef(knon, Kcoefhq(:,:), gamaq(:,:), delp(:,:), q(:,:), &
         Ccoef_Q(:,:), Dcoef_Q(:,:), Acoef_Q, Bcoef_Q)

!****************************************************************************************
! 5)
! Calculte the coefficients C and D for potentiel entalpie, H
!
!****************************************************************************************
    h_old(:,:) = 0.0

    DO k=1,klev
       DO i = 1, knon
          ! convertie la temperature en entalpie potentielle
          h_old(i,k) = RCPD * temp(i,k) * &
               (psref(i)/pplay(i,k))**RKAPPA
       ENDDO
    ENDDO

    CALL calc_coef(knon, Kcoefhq(:,:), gamah(:,:), delp(:,:), h_old(:,:), &
         Ccoef_H(:,:), Dcoef_H(:,:), Acoef_H, Bcoef_H)

!****************************************************************************************
! 6)
! Return the first layer in output variables
!
!****************************************************************************************
    Acoef_H_out = Acoef_H
    Bcoef_H_out = Bcoef_H
    Acoef_Q_out = Acoef_Q
    Bcoef_Q_out = Bcoef_Q

!****************************************************************************************
! 7)
! If Pbl is split, return also the other layers in output variables
!
!****************************************************************************************
!!! jyg le 07/02/2012
!!jyg       IF (mod(iflag_pbl_split,2) .eq.1) THEN
       IF (mod(iflag_pbl_split,10) .ge.1) THEN
!!! nrlmd le 02/05/2011
    DO k= 1, klev
      DO i= 1, klon
        Ccoef_H_out(i,k) = Ccoef_H(i,k)
        Dcoef_H_out(i,k) = Dcoef_H(i,k)
        Ccoef_Q_out(i,k) = Ccoef_Q(i,k)
        Dcoef_Q_out(i,k) = Dcoef_Q(i,k)
        Kcoef_hq_out(i,k) = Kcoefhq(i,k)
          IF (k.eq.1) THEN
            gama_h_out(i,k)  = 0.
            gama_q_out(i,k)  = 0.
          ELSE
            gama_h_out(i,k)  = gamah(i,k)
            gama_q_out(i,k)  = gamaq(i,k)
          ENDIF
      ENDDO
    ENDDO
!!!
       ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
!!!

  END SUBROUTINE climb_hq_down
!
!****************************************************************************************
!
  SUBROUTINE calc_coef(knon, Kcoef, gama, delp, X, Ccoef, Dcoef, Acoef, Bcoef)
!
! Calculate the coefficients C and D in : X(k) = C(k) + D(k)*X(k-1)
! where X is H or Q, and k the vertical level k=1,klev
!
    INCLUDE "YOMCST.h"
! Input arguments
!****************************************************************************************
    INTEGER, INTENT(IN)                      :: knon
    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Kcoef, delp
    REAL, DIMENSION(klon,klev), INTENT(IN)   :: X
    REAL, DIMENSION(klon,2:klev), INTENT(IN) :: gama

! Output arguments
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(OUT)       :: Acoef, Bcoef
    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: Ccoef, Dcoef

! Local variables
!****************************************************************************************
    INTEGER                                  :: k, i
    REAL                                     :: buf

!****************************************************************************************
! Niveau au sommet, k=klev
!
!****************************************************************************************
    Ccoef(:,:) = 0.0
    Dcoef(:,:) = 0.0

    DO i = 1, knon
       buf = delp(i,klev) + Kcoef(i,klev)

       Ccoef(i,klev) = (X(i,klev)*delp(i,klev) - Kcoef(i,klev)*gama(i,klev))/buf
       Dcoef(i,klev) = Kcoef(i,klev)/buf
    END DO


!****************************************************************************************
! Niveau  (klev-1) <= k <= 2
!
!****************************************************************************************

    DO k=(klev-1),2,-1
       DO i = 1, knon
          buf = delp(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1))
          Ccoef(i,k) = (X(i,k)*delp(i,k) + Kcoef(i,k+1)*Ccoef(i,k+1) + &
               Kcoef(i,k+1)*gama(i,k+1) - Kcoef(i,k)*gama(i,k))/buf
          Dcoef(i,k) = Kcoef(i,k)/buf
       END DO
    END DO

!****************************************************************************************
! Niveau k=1
!
!****************************************************************************************

    DO i = 1, knon
       buf = delp(i,1) + Kcoef(i,2)*(1.-Dcoef(i,2))
       Acoef(i) = (X(i,1)*delp(i,1) + Kcoef(i,2)*(gama(i,2)+Ccoef(i,2)))/buf
       Bcoef(i) = -1. * RG / buf
    END DO

  END SUBROUTINE calc_coef
!
!****************************************************************************************
!
  SUBROUTINE climb_hq_up(knon, dtime, t_old, q_old, &
       flx_q1, flx_h1, paprs, pplay, &
!!! nrlmd le 02/05/2011
       Acoef_H_in, Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in, &
       Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in, &
       Kcoef_hq_in, gama_q_in, gama_h_in, &
!!!
       flux_q, flux_h, d_q, d_t)
!
! This routine calculates the flux and tendency of the specific humidity q and
! the potential engergi H.
! The quantities q and H are calculated according to
! X(k) = C(k) + D(k)*X(k-1) for X=[q,H], where the coefficients
! C and D are known from before and k is index of the vertical layer.
!

! Input arguments
!****************************************************************************************
    INTEGER, INTENT(IN)                      :: knon
    REAL, INTENT(IN)                         :: dtime
    REAL, DIMENSION(klon,klev), INTENT(IN)   :: t_old, q_old
    REAL, DIMENSION(klon), INTENT(IN)        :: flx_q1, flx_h1
    REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
    REAL, DIMENSION(klon,klev), INTENT(IN)   :: pplay

!!! nrlmd le 02/05/2011
    REAL, DIMENSION(klon), INTENT(IN)        :: Acoef_H_in,Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in
    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in
    REAL, DIMENSION(klon,klev), INTENT(IN)   :: Kcoef_hq_in, gama_q_in, gama_h_in
!!!

! Output arguments
!****************************************************************************************
    REAL, DIMENSION(klon,klev), INTENT(OUT)  :: flux_q, flux_h, d_q, d_t

! Local variables
!****************************************************************************************
    LOGICAL, SAVE                            :: last=.FALSE.
    REAL, DIMENSION(klon,klev)               :: h_new, q_new
    REAL, DIMENSION(klon)                    :: psref
    INTEGER                                  :: k, i, ierr

! Include
!****************************************************************************************
    INCLUDE "YOMCST.h"
    INCLUDE "compbl.h"

!****************************************************************************************
! 1)
! Definition of some variables
    REAL, DIMENSION(klon,klev)               :: d_h, zairm
!
!****************************************************************************************
    flux_q(:,:) = 0.0
    flux_h(:,:) = 0.0
    d_q(:,:)    = 0.0
    d_t(:,:)    = 0.0
    d_h(:,:)    = 0.0
    f_h_bnd(:)= 0.0

    psref(1:knon) = paprs(1:knon,1)

!!! jyg le 07/02/2012
!!jyg       IF (mod(iflag_pbl_split,2) .eq.1) THEN
       IF (mod(iflag_pbl_split,10) .ge.1) THEN
!!! nrlmd le 02/05/2011
    DO i = 1, knon
      Acoef_H(i)=Acoef_H_in(i)
      Acoef_Q(i)=Acoef_Q_in(i)
      Bcoef_H(i)=Bcoef_H_in(i)
      Bcoef_Q(i)=Bcoef_Q_in(i)
    ENDDO
    DO k = 1, klev
      DO i = 1, knon
        Ccoef_H(i,k)=Ccoef_H_in(i,k)
        Ccoef_Q(i,k)=Ccoef_Q_in(i,k)
        Dcoef_H(i,k)=Dcoef_H_in(i,k)
        Dcoef_Q(i,k)=Dcoef_Q_in(i,k)
        Kcoefhq(i,k)=Kcoef_hq_in(i,k)
          IF (k.gt.1) THEN
            gamah(i,k)=gama_h_in(i,k)
            gamaq(i,k)=gama_q_in(i,k)
          ENDIF
      ENDDO
    ENDDO
!!!
       ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
!!!

!****************************************************************************************
! 2)
! Calculation of Q and H
!
!****************************************************************************************

!- First layer
    q_new(1:knon,1) = Acoef_Q(1:knon) + Bcoef_Q(1:knon)*flx_q1(1:knon)*dtime
    h_new(1:knon,1) = Acoef_H(1:knon) + Bcoef_H(1:knon)*flx_h1(1:knon)*dtime
    f_h_bnd(1:knon) = flx_h1(1:knon)
!- All the other layers
    DO k = 2, klev
       DO i = 1, knon
          q_new(i,k) = Ccoef_Q(i,k) + Dcoef_Q(i,k)*q_new(i,k-1)
          h_new(i,k) = Ccoef_H(i,k) + Dcoef_H(i,k)*h_new(i,k-1)
       END DO
    END DO
!****************************************************************************************
! 3)
! Calculation of the flux for Q and H
!
!****************************************************************************************

!- The flux at first layer, k=1
    flux_q(1:knon,1)=flx_q1(1:knon)
    flux_h(1:knon,1)=flx_h1(1:knon)

!- The flux at all layers above surface
    DO k = 2, klev
       DO i = 1, knon
          flux_q(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
               (q_new(i,k)-q_new(i,k-1)+gamaq(i,k))

          flux_h(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
               (h_new(i,k)-h_new(i,k-1)+gamah(i,k))
       END DO
    END DO

!****************************************************************************************
! 4)
! Calculation of tendency for Q and H
!
!****************************************************************************************
    d_h_col_vdf(:) = 0.0
    DO k = 1, klev
       DO i = 1, knon
          d_t(i,k) = h_new(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD - t_old(i,k)
          d_q(i,k) = q_new(i,k) - q_old(i,k)
          d_h(i,k) = h_new(i,k) - h_old(i,k)
!JLD          d_t(i,k) = d_h(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD    !correction a venir
    ! layer air mass
          zairm(i, k) = (paprs(i,k)-paprs(i,k+1))/rg
          d_h_col_vdf(i) = d_h_col_vdf(i) + d_h(i,k)*zairm(i,k)
        END DO
    END DO

!****************************************************************************************
! Some deallocations
!
!****************************************************************************************
    IF (last) THEN
       DEALLOCATE(Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H,stat=ierr)
       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H, ierr=', ierr
       DEALLOCATE(Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H,stat=ierr)
       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H, ierr=', ierr
       DEALLOCATE(gamaq, gamah,stat=ierr)
       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate gamaq, gamah, ierr=', ierr
       DEALLOCATE(Kcoefhq,stat=ierr)
       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Kcoefhq, ierr=', ierr
       DEALLOCATE(h_old, d_h_col_vdf, f_h_bnd, stat=ierr)
       IF ( ierr /= 0 )  PRINT*,' pb in dealllocate h_old, d_h_col_vdf, f_h_bnd, ierr=', ierr
    END IF
  END SUBROUTINE climb_hq_up
!
!****************************************************************************************
!
END MODULE climb_hq_mod








Generated by: GCOVR (Version 4.2)

Line	Branch	Exec	Source
1			MODULE climb_hq_mod
2			!
3			! Module to solve the verctical diffusion of "q" and "H";
4			! specific humidity and potential energi.
5			!
6			USE dimphy
7
8			IMPLICIT NONE
9			SAVE
10			PRIVATE
11			PUBLIC :: climb_hq_down, climb_hq_up, d_h_col_vdf, f_h_bnd
12
13			REAL, DIMENSION(:,:), ALLOCATABLE :: gamaq, gamah
14			!$OMP THREADPRIVATE(gamaq,gamah)
15			REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_Q, Dcoef_Q
16			!$OMP THREADPRIVATE(Ccoef_Q, Dcoef_Q)
17			REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_H, Dcoef_H
18			!$OMP THREADPRIVATE(Ccoef_H, Dcoef_H)
19			REAL, DIMENSION(:), ALLOCATABLE :: Acoef_Q, Bcoef_Q
20			!$OMP THREADPRIVATE(Acoef_Q, Bcoef_Q)
21			REAL, DIMENSION(:), ALLOCATABLE :: Acoef_H, Bcoef_H
22			!$OMP THREADPRIVATE(Acoef_H, Bcoef_H)
23			REAL, DIMENSION(:,:), ALLOCATABLE :: Kcoefhq
24			!$OMP THREADPRIVATE(Kcoefhq)
25			REAL, SAVE, DIMENSION(:,:), ALLOCATABLE :: h_old ! for diagnostics, h before solving diffusion
26			!$OMP THREADPRIVATE(h_old)
27			REAL, SAVE, DIMENSION(:), ALLOCATABLE :: d_h_col_vdf ! for diagnostics, vertical integral of enthalpy change
28			!$OMP THREADPRIVATE(d_h_col_vdf)
29			REAL, SAVE, DIMENSION(:), ALLOCATABLE :: f_h_bnd ! for diagnostics, enthalpy flux at surface
30			!$OMP THREADPRIVATE(f_h_bnd)
31
32			CONTAINS
33			!
34			!****************************************************************************************
35			!
36		36408058	SUBROUTINE climb_hq_down(knon, coefhq, paprs, pplay, &
37			delp, temp, q, dtime, &
38			!!! nrlmd le 02/05/2011
39			Ccoef_H_out, Ccoef_Q_out, Dcoef_H_out, Dcoef_Q_out, &
40			Kcoef_hq_out, gama_q_out, gama_h_out, &
41			!!!
42		1152	Acoef_H_out, Acoef_Q_out, Bcoef_H_out, Bcoef_Q_out)
43
44			! This routine calculates recursivly the coefficients C and D
45			! for the quantity X=[Q,H] in equation X(k) = C(k) + D(k)*X(k-1), where k is
46			! the index of the vertical layer.
47			!
48			! Input arguments
49			!****************************************************************************************
50			INTEGER, INTENT(IN) :: knon
51			REAL, DIMENSION(klon,klev), INTENT(IN) :: coefhq
52			REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay
53			REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
54			REAL, DIMENSION(klon,klev), INTENT(IN) :: temp, delp ! temperature
55			REAL, DIMENSION(klon,klev), INTENT(IN) :: q
56			REAL, INTENT(IN) :: dtime
57
58			! Output arguments
59			!****************************************************************************************
60			REAL, DIMENSION(klon), INTENT(OUT) :: Acoef_H_out
61			REAL, DIMENSION(klon), INTENT(OUT) :: Acoef_Q_out
62			REAL, DIMENSION(klon), INTENT(OUT) :: Bcoef_H_out
63			REAL, DIMENSION(klon), INTENT(OUT) :: Bcoef_Q_out
64
65			!!! nrlmd le 02/05/2011
66			REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef_H_out
67			REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef_Q_out
68			REAL, DIMENSION(klon,klev), INTENT(OUT) :: Dcoef_H_out
69			REAL, DIMENSION(klon,klev), INTENT(OUT) :: Dcoef_Q_out
70			REAL, DIMENSION(klon,klev), INTENT(OUT) :: Kcoef_hq_out
71			REAL, DIMENSION(klon,klev), INTENT(OUT) :: gama_q_out
72			REAL, DIMENSION(klon,klev), INTENT(OUT) :: gama_h_out
73			!!!
74
75			! Local variables
76			!****************************************************************************************
77			LOGICAL, SAVE :: first=.TRUE.
78			!$OMP THREADPRIVATE(first)
79			! JLD now renamed h_old and declared in module
80			! REAL, DIMENSION(klon,klev) :: local_H
81		2304	REAL, DIMENSION(klon) :: psref
82			REAL :: delz, pkh
83			INTEGER :: k, i, ierr
84			! Include
85			!****************************************************************************************
86			INCLUDE "YOMCST.h"
87			INCLUDE "compbl.h"
88
89
90			!****************************************************************************************
91			! 1)
92			! Allocation at first time step only
93			!
94			!****************************************************************************************
95
96	✓✓	1152	IF (first) THEN
97		1	first=.FALSE.
98	✓✗✓✗ ✓✗✓✗ ✗✓✓✗	2	ALLOCATE(Ccoef_Q(klon,klev), STAT=ierr)
99	✗✓	1	IF ( ierr /= 0 ) PRINT*,' pb in allloc Ccoef_Q, ierr=', ierr
100
101	✓✗✓✗ ✓✗✓✗ ✗✓✓✗	2	ALLOCATE(Dcoef_Q(klon,klev), STAT=ierr)
102	✗✓	1	IF ( ierr /= 0 ) PRINT*,' pb in allloc Dcoef_Q, ierr=', ierr
103
104	✓✗✓✗ ✓✗✓✗ ✗✓✓✗	2	ALLOCATE(Ccoef_H(klon,klev), STAT=ierr)
105	✗✓	1	IF ( ierr /= 0 ) PRINT*,' pb in allloc Ccoef_H, ierr=', ierr
106
107	✓✗✓✗ ✓✗✓✗ ✗✓✓✗	2	ALLOCATE(Dcoef_H(klon,klev), STAT=ierr)
108	✗✓	1	IF ( ierr /= 0 ) PRINT*,' pb in allloc Dcoef_H, ierr=', ierr
109
110	✓✗✓✗ ✗✓✓✗ ✓✗✓✗ ✗✓✓✗ ✓✗✓✗ ✗✓✓✗ ✓✗✓✗ ✗✓✓✗	1	ALLOCATE(Acoef_Q(klon), Bcoef_Q(klon), Acoef_H(klon), Bcoef_H(klon), STAT=ierr)
111	✗✓	1	IF ( ierr /= 0 ) PRINT*,' pb in allloc Acoef_X and Bcoef_X, ierr=', ierr
112
113	✓✗✓✗ ✓✗✓✗ ✗✓✓✗	2	ALLOCATE(Kcoefhq(klon,klev), STAT=ierr)
114	✗✓	1	IF ( ierr /= 0 ) PRINT*,' pb in allloc Kcoefhq, ierr=', ierr
115
116	✓✗✓✗ ✓✗✓✗ ✗✓✓✗	2	ALLOCATE(gamaq(1:klon,2:klev), STAT=ierr)
117	✗✓	1	IF ( ierr /= 0 ) PRINT*,' pb in allloc gamaq, ierr=', ierr
118
119	✓✗✓✗ ✓✗✓✗ ✗✓✓✗	2	ALLOCATE(gamah(1:klon,2:klev), STAT=ierr)
120	✗✓	1	IF ( ierr /= 0 ) PRINT*,' pb in allloc gamah, ierr=', ierr
121
122	✓✗✓✗ ✓✗✓✗ ✗✓✓✗	2	ALLOCATE(h_old(klon,klev), STAT=ierr)
123	✗✓	1	IF ( ierr /= 0 ) PRINT*,' pb in allloc h_old, ierr=', ierr
124
125	✓✗✓✗ ✗✓✓✗	1	ALLOCATE(d_h_col_vdf(klon), STAT=ierr)
126	✗✓	1	IF ( ierr /= 0 ) PRINT*,' pb in allloc d_h_col_vdf, ierr=', ierr
127
128	✓✗✓✗ ✗✓✓✗	1	ALLOCATE(f_h_bnd(klon), STAT=ierr)
129	✗✓	1	IF ( ierr /= 0 ) PRINT*,' pb in allloc f_h_bnd, ierr=', ierr
130			END IF
131
132			!****************************************************************************************
133			! 2)
134			! Definition of the coeficient K
135			!
136			!****************************************************************************************
137	✓✓✓✓	44704512	Kcoefhq(:,:) = 0.0
138	✓✓	44928	DO k = 2, klev
139	✓✓	18011404	DO i = 1, knon
140			Kcoefhq(i,k) = &
141			coefhq(i,k)RGRG*dtime /(pplay(i,k-1)-pplay(i,k)) &
142		43776	(paprs(i,k)2/(temp(i,k)+temp(i,k-1))/RD)**2
143			ENDDO
144			ENDDO
145
146			!****************************************************************************************
147			! 3)
148			! Calculation of gama for "Q" and "H"
149			!
150			!****************************************************************************************
151			! surface pressure is used as reference
152	✓✓	1146240	psref(:) = paprs(:,1)
153
154			! definition of gama
155		1152	IF (iflag_pbl == 1) THEN
156			gamaq(:,:) = 0.0
157			gamah(:,:) = -1.0e-03
158			gamah(:,2) = -2.5e-03
159
160			! conversion de gama
161			DO k = 2, klev
162			DO i = 1, knon
163			delz = RD * (temp(i,k-1)+temp(i,k)) / &
164			2.0 / RG / paprs(i,k) * (pplay(i,k-1)-pplay(i,k))
165			pkh = (psref(i)/paprs(i,k))**RKAPPA
166
167			! convertie gradient verticale d'humidite specifique en difference d'humidite specifique entre centre de couches
168			gamaq(i,k) = gamaq(i,k) * delz
169			! convertie gradient verticale de temperature en difference de temperature potentielle entre centre de couches
170			gamah(i,k) = gamah(i,k) * delz * RCPD * pkh
171			ENDDO
172			ENDDO
173
174			ELSE
175	✓✓✓✓	43558272	gamaq(:,:) = 0.0
176	✓✓✓✓	43558272	gamah(:,:) = 0.0
177			ENDIF
178
179
180			!****************************************************************************************
181			! 4)
182			! Calculte the coefficients C and D for specific humidity, q
183			!
184			!****************************************************************************************
185
186			CALL calc_coef(knon, Kcoefhq(:,:), gamaq(:,:), delp(:,:), q(:,:), &
187		1152	Ccoef_Q(:,:), Dcoef_Q(:,:), Acoef_Q, Bcoef_Q)
188
189			!****************************************************************************************
190			! 5)
191			! Calculte the coefficients C and D for potentiel entalpie, H
192			!
193			!****************************************************************************************
194	✓✓✓✓	44704512	h_old(:,:) = 0.0
195
196	✓✓	46080	DO k=1,klev
197	✓✓	18485358	DO i = 1, knon
198			! convertie la temperature en entalpie potentielle
199			h_old(i,k) = RCPD * temp(i,k) * &
200		44928	(psref(i)/pplay(i,k))**RKAPPA
201			ENDDO
202			ENDDO
203
204			CALL calc_coef(knon, Kcoefhq(:,:), gamah(:,:), delp(:,:), h_old(:,:), &
205		1152	Ccoef_H(:,:), Dcoef_H(:,:), Acoef_H, Bcoef_H)
206
207			!****************************************************************************************
208			! 6)
209			! Return the first layer in output variables
210			!
211			!****************************************************************************************
212	✓✓	1146240	Acoef_H_out = Acoef_H
213	✓✓	1146240	Bcoef_H_out = Bcoef_H
214	✓✓	1146240	Acoef_Q_out = Acoef_Q
215	✓✓	1146240	Bcoef_Q_out = Bcoef_Q
216
217			!****************************************************************************************
218			! 7)
219			! If Pbl is split, return also the other layers in output variables
220			!
221			!****************************************************************************************
222			!!! jyg le 07/02/2012
223			!!jyg IF (mod(iflag_pbl_split,2) .eq.1) THEN
224		1152	IF (mod(iflag_pbl_split,10) .ge.1) THEN
225			!!! nrlmd le 02/05/2011
226			DO k= 1, klev
227			DO i= 1, klon
228			Ccoef_H_out(i,k) = Ccoef_H(i,k)
229			Dcoef_H_out(i,k) = Dcoef_H(i,k)
230			Ccoef_Q_out(i,k) = Ccoef_Q(i,k)
231			Dcoef_Q_out(i,k) = Dcoef_Q(i,k)
232			Kcoef_hq_out(i,k) = Kcoefhq(i,k)
233			IF (k.eq.1) THEN
234			gama_h_out(i,k) = 0.
235			gama_q_out(i,k) = 0.
236			ELSE
237			gama_h_out(i,k) = gamah(i,k)
238			gama_q_out(i,k) = gamaq(i,k)
239			ENDIF
240			ENDDO
241			ENDDO
242			!!!
243			ENDIF ! (mod(iflag_pbl_split,2) .ge.1)
244			!!!
245
246		1152	END SUBROUTINE climb_hq_down
247			!
248			!****************************************************************************************
249			!
250		2304	SUBROUTINE calc_coef(knon, Kcoef, gama, delp, X, Ccoef, Dcoef, Acoef, Bcoef)
251			!
252			! Calculate the coefficients C and D in : X(k) = C(k) + D(k)*X(k-1)
253			! where X is H or Q, and k the vertical level k=1,klev
254			!
255			INCLUDE "YOMCST.h"
256			! Input arguments
257			!****************************************************************************************
258			INTEGER, INTENT(IN) :: knon
259			REAL, DIMENSION(klon,klev), INTENT(IN) :: Kcoef, delp
260			REAL, DIMENSION(klon,klev), INTENT(IN) :: X
261			REAL, DIMENSION(klon,2:klev), INTENT(IN) :: gama
262
263			! Output arguments
264			!****************************************************************************************
265			REAL, DIMENSION(klon), INTENT(OUT) :: Acoef, Bcoef
266			REAL, DIMENSION(klon,klev), INTENT(OUT) :: Ccoef, Dcoef
267
268			! Local variables
269			!****************************************************************************************
270			INTEGER :: k, i
271			REAL :: buf
272
273			!****************************************************************************************
274			! Niveau au sommet, k=klev
275			!
276			!****************************************************************************************
277	✓✓✓✓	89409024	Ccoef(:,:) = 0.0
278	✓✓✓✓	89409024	Dcoef(:,:) = 0.0
279
280	✓✓	947908	DO i = 1, knon
281		945604	buf = delp(i,klev) + Kcoef(i,klev)
282
283		945604	Ccoef(i,klev) = (X(i,klev)delp(i,klev) - Kcoef(i,klev)gama(i,klev))/buf
284		947908	Dcoef(i,klev) = Kcoef(i,klev)/buf
285			END DO
286
287
288			!****************************************************************************************
289			! Niveau (klev-1) <= k <= 2
290			!
291			!****************************************************************************************
292
293	✓✓	87552	DO k=(klev-1),2,-1
294	✓✓	35074900	DO i = 1, knon
295		34987348	buf = delp(i,k) + Kcoef(i,k) + Kcoef(i,k+1)*(1.-Dcoef(i,k+1))
296			Ccoef(i,k) = (X(i,k)delp(i,k) + Kcoef(i,k+1)Ccoef(i,k+1) + &
297		34987348	Kcoef(i,k+1)gama(i,k+1) - Kcoef(i,k)gama(i,k))/buf
298		35072596	Dcoef(i,k) = Kcoef(i,k)/buf
299			END DO
300			END DO
301
302			!****************************************************************************************
303			! Niveau k=1
304			!
305			!****************************************************************************************
306
307	✓✓	947908	DO i = 1, knon
308		945604	buf = delp(i,1) + Kcoef(i,2)*(1.-Dcoef(i,2))
309		945604	Acoef(i) = (X(i,1)delp(i,1) + Kcoef(i,2)(gama(i,2)+Ccoef(i,2)))/buf
310		2304	Bcoef(i) = -1. * RG / buf
311			END DO
312
313	✗✓✗✓	36408058	END SUBROUTINE calc_coef
314			!
315			!****************************************************************************************
316			!
317		1152	SUBROUTINE climb_hq_up(knon, dtime, t_old, q_old, &
318			flx_q1, flx_h1, paprs, pplay, &
319			!!! nrlmd le 02/05/2011
320		1152	Acoef_H_in, Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in, &
321			Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in, &
322			Kcoef_hq_in, gama_q_in, gama_h_in, &
323			!!!
324		1152	flux_q, flux_h, d_q, d_t)
325			!
326			! This routine calculates the flux and tendency of the specific humidity q and
327			! the potential engergi H.
328			! The quantities q and H are calculated according to
329			! X(k) = C(k) + D(k)*X(k-1) for X=[q,H], where the coefficients
330			! C and D are known from before and k is index of the vertical layer.
331			!
332
333			! Input arguments
334			!****************************************************************************************
335			INTEGER, INTENT(IN) :: knon
336			REAL, INTENT(IN) :: dtime
337			REAL, DIMENSION(klon,klev), INTENT(IN) :: t_old, q_old
338			REAL, DIMENSION(klon), INTENT(IN) :: flx_q1, flx_h1
339			REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs
340			REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay
341
342			!!! nrlmd le 02/05/2011
343			REAL, DIMENSION(klon), INTENT(IN) :: Acoef_H_in,Acoef_Q_in, Bcoef_H_in, Bcoef_Q_in
344			REAL, DIMENSION(klon,klev), INTENT(IN) :: Ccoef_H_in, Ccoef_Q_in, Dcoef_H_in, Dcoef_Q_in
345			REAL, DIMENSION(klon,klev), INTENT(IN) :: Kcoef_hq_in, gama_q_in, gama_h_in
346			!!!
347
348			! Output arguments
349			!****************************************************************************************
350			REAL, DIMENSION(klon,klev), INTENT(OUT) :: flux_q, flux_h, d_q, d_t
351
352			! Local variables
353			!****************************************************************************************
354			LOGICAL, SAVE :: last=.FALSE.
355		2304	REAL, DIMENSION(klon,klev) :: h_new, q_new
356		2304	REAL, DIMENSION(klon) :: psref
357			INTEGER :: k, i, ierr
358
359			! Include
360			!****************************************************************************************
361			INCLUDE "YOMCST.h"
362			INCLUDE "compbl.h"
363
364			!****************************************************************************************
365			! 1)
366			! Definition of some variables
367		1152	REAL, DIMENSION(klon,klev) :: d_h, zairm
368			!
369			!****************************************************************************************
370	✓✓✓✓	44704512	flux_q(:,:) = 0.0
371	✓✓✓✓	44704512	flux_h(:,:) = 0.0
372	✓✓✓✓	44704512	d_q(:,:) = 0.0
373	✓✓✓✓	44704512	d_t(:,:) = 0.0
374	✓✓✓✓	44704512	d_h(:,:) = 0.0
375	✓✓	1146240	f_h_bnd(:)= 0.0
376
377	✓✓	473954	psref(1:knon) = paprs(1:knon,1)
378
379			!!! jyg le 07/02/2012
380			!!jyg IF (mod(iflag_pbl_split,2) .eq.1) THEN
381	✗✓	1152	IF (mod(iflag_pbl_split,10) .ge.1) THEN
382			!!! nrlmd le 02/05/2011
383			DO i = 1, knon
384			Acoef_H(i)=Acoef_H_in(i)
385			Acoef_Q(i)=Acoef_Q_in(i)
386			Bcoef_H(i)=Bcoef_H_in(i)
387			Bcoef_Q(i)=Bcoef_Q_in(i)
388			ENDDO
389			DO k = 1, klev
390			DO i = 1, knon
391			Ccoef_H(i,k)=Ccoef_H_in(i,k)
392			Ccoef_Q(i,k)=Ccoef_Q_in(i,k)
393			Dcoef_H(i,k)=Dcoef_H_in(i,k)
394			Dcoef_Q(i,k)=Dcoef_Q_in(i,k)
395			Kcoefhq(i,k)=Kcoef_hq_in(i,k)
396			IF (k.gt.1) THEN
397			gamah(i,k)=gama_h_in(i,k)
398			gamaq(i,k)=gama_q_in(i,k)
399			ENDIF
400			ENDDO
401			ENDDO
402			!!!
403			ENDIF ! (mod(iflag_pbl_split,2) .ge.1)
404			!!!
405
406			!****************************************************************************************
407			! 2)
408			! Calculation of Q and H
409			!
410			!****************************************************************************************
411
412			!- First layer
413	✓✓	473954	q_new(1:knon,1) = Acoef_Q(1:knon) + Bcoef_Q(1:knon)flx_q1(1:knon)dtime
414	✓✓	473954	h_new(1:knon,1) = Acoef_H(1:knon) + Bcoef_H(1:knon)flx_h1(1:knon)dtime
415	✓✓	473954	f_h_bnd(1:knon) = flx_h1(1:knon)
416			!- All the other layers
417	✓✓	44928	DO k = 2, klev
418	✓✓	18011404	DO i = 1, knon
419		17966476	q_new(i,k) = Ccoef_Q(i,k) + Dcoef_Q(i,k)*q_new(i,k-1)
420		18010252	h_new(i,k) = Ccoef_H(i,k) + Dcoef_H(i,k)*h_new(i,k-1)
421			END DO
422			END DO
423			!****************************************************************************************
424			! 3)
425			! Calculation of the flux for Q and H
426			!
427			!****************************************************************************************
428
429			!- The flux at first layer, k=1
430	✓✓	473954	flux_q(1:knon,1)=flx_q1(1:knon)
431	✓✓	473954	flux_h(1:knon,1)=flx_h1(1:knon)
432
433			!- The flux at all layers above surface
434	✓✓	44928	DO k = 2, klev
435	✓✓	18011404	DO i = 1, knon
436			flux_q(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
437		17966476	(q_new(i,k)-q_new(i,k-1)+gamaq(i,k))
438
439			flux_h(i,k) = (Kcoefhq(i,k)/RG/dtime) * &
440		18010252	(h_new(i,k)-h_new(i,k-1)+gamah(i,k))
441			END DO
442			END DO
443
444			!****************************************************************************************
445			! 4)
446			! Calculation of tendency for Q and H
447			!
448			!****************************************************************************************
449	✓✓	1146240	d_h_col_vdf(:) = 0.0
450	✓✓	46080	DO k = 1, klev
451	✓✓	18485358	DO i = 1, knon
452		18439278	d_t(i,k) = h_new(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD - t_old(i,k)
453		18439278	d_q(i,k) = q_new(i,k) - q_old(i,k)
454		18439278	d_h(i,k) = h_new(i,k) - h_old(i,k)
455			!JLD d_t(i,k) = d_h(i,k)/(psref(i)/pplay(i,k))**RKAPPA/RCPD !correction a venir
456			! layer air mass
457		18439278	zairm(i, k) = (paprs(i,k)-paprs(i,k+1))/rg
458		18484206	d_h_col_vdf(i) = d_h_col_vdf(i) + d_h(i,k)*zairm(i,k)
459			END DO
460			END DO
461
462			!****************************************************************************************
463			! Some deallocations
464			!
465			!****************************************************************************************
466	✗✓	1152	IF (last) THEN
467			DEALLOCATE(Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H,stat=ierr)
468			IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Ccoef_Q, Dcoef_Q, Ccoef_H, Dcoef_H, ierr=', ierr
469			DEALLOCATE(Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H,stat=ierr)
470			IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Acoef_Q, Bcoef_Q, Acoef_H, Bcoef_H, ierr=', ierr
471			DEALLOCATE(gamaq, gamah,stat=ierr)
472			IF ( ierr /= 0 ) PRINT*,' pb in dealllocate gamaq, gamah, ierr=', ierr
473			DEALLOCATE(Kcoefhq,stat=ierr)
474			IF ( ierr /= 0 ) PRINT*,' pb in dealllocate Kcoefhq, ierr=', ierr
475			DEALLOCATE(h_old, d_h_col_vdf, f_h_bnd, stat=ierr)
476			IF ( ierr /= 0 ) PRINT*,' pb in dealllocate h_old, d_h_col_vdf, f_h_bnd, ierr=', ierr
477			END IF
478		945604	END SUBROUTINE climb_hq_up
479			!
480			!****************************************************************************************
481			!
482			END MODULE climb_hq_mod
483
484
485
486
487
488