climb_hq_mod.F90

Go to the documentation of this file.

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

  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)

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)
    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
    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 
!
!****************************************************************************************
    local_h(:,:) = 0.0

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

    CALL calc_coef(knon, kcoefhq(:,:), gamah(:,:), delp(:,:), local_h(:,:), &
         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
       IF (mod(iflag_pbl_split,2) .eq.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) .eq.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
!
!****************************************************************************************
    flux_q(:,:) = 0.0
    flux_h(:,:) = 0.0
    d_q(:,:)    = 0.0
    d_t(:,:)    = 0.0

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

!!! jyg le 07/02/2012
       IF (mod(iflag_pbl_split,2) .eq.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) .eq.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
    
!- 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
!
!****************************************************************************************

    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)
       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
    END IF
  END SUBROUTINE climb_hq_up
!
!****************************************************************************************
!
END MODULE climb_hq_mod