MODULE climb_qbs_mod
!
! Module to solve the verctical diffusion of blowing snow; 
!

  IMPLICIT NONE
  SAVE 
  PRIVATE
  PUBLIC :: climb_qbs_down, climb_qbs_up, climb_qbs_init, climb_qbs_finalize

  REAL, DIMENSION(:,:), ALLOCATABLE :: gamaqbs
  !$OMP THREADPRIVATE(gamaqbs)
  REAL, DIMENSION(:,:), ALLOCATABLE :: Ccoef_QBS, Dcoef_QBS
  !$OMP THREADPRIVATE(Ccoef_QBS, Dcoef_QBS)
  REAL, DIMENSION(:), ALLOCATABLE   :: Acoef_QBS, Bcoef_QBS
  !$OMP THREADPRIVATE(Acoef_QBS, Bcoef_QBS)
  REAL, DIMENSION(:,:), ALLOCATABLE :: Kcoefqbs
  !$OMP THREADPRIVATE(Kcoefqbs)

  LOGICAL, SAVE                            :: first=.TRUE.
  !$OMP THREADPRIVATE(first)
 
CONTAINS
  SUBROUTINE climb_qbs_init
  USE dimphy, ONLY : klon, klev
  IMPLICIT NONE 
    INTEGER :: ierr
!****************************************************************************************
! 1)
! Allocation at first time step only
!   
!****************************************************************************************

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

       ALLOCATE(Acoef_QBS(klon), Bcoef_QBS(klon), STAT=ierr)
       IF ( ierr /= 0 )  PRINT*,' pb in allloc Acoef_BS and Bcoef_BS, ierr=', ierr
       Acoef_QBS(:) = 0. ; Bcoef_QBS(:) = 0.

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

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

    END IF

  END SUBROUTINE climb_qbs_init

  SUBROUTINE climb_qbs_finalize
  IMPLICIT NONE
  INTEGER :: ierr
    !****************************************************************************************
    ! Some deallocations
    !
    !****************************************************************************************
    DEALLOCATE(Ccoef_QBS, Dcoef_QBS,stat=ierr)    
    IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Ccoef_QBS, Dcoef_QBS, ierr=', ierr
    DEALLOCATE(Acoef_QBS, Bcoef_QBS,stat=ierr)    
    IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Acoef_QBS, Bcoef_QBS, ierr=', ierr
    DEALLOCATE(gamaqbs,stat=ierr)
    IF ( ierr /= 0 )  PRINT*,' pb in dealllocate gamaqbs, ierr=', ierr
    DEALLOCATE(Kcoefqbs,stat=ierr)
    IF ( ierr /= 0 )  PRINT*,' pb in dealllocate Kcoefqbs, ierr=', ierr

  END SUBROUTINE climb_qbs_finalize

!
!****************************************************************************************
!
  SUBROUTINE climb_qbs_down(knon, ni, coefqbs, paprs, pplay, &
       delp, temp, qbs, dtime, &
       Ccoef_QBS_out, Dcoef_QBS_out, &
       Kcoef_qbs_out, gama_qbs_out, &
       Acoef_QBS_out, Bcoef_QBS_out)
!$gpum horizontal knon

! This routine calculates recursivly the coefficients C and D
! for the quantity X=[QBS] in equation X(k) = C(k) + D(k)*X(k-1), where k is
! the index of the vertical layer.
USE yomcst_mod_h
USE compbl_mod_h
USE dimphy, ONLY : klev
! Input arguments
!****************************************************************************************
    INTEGER, INTENT(IN)                      :: knon
    INTEGER, INTENT(IN)                      :: ni(knon)
    REAL, DIMENSION(knon,klev), INTENT(IN)   :: coefqbs
    REAL, DIMENSION(knon,klev), INTENT(IN)   :: pplay 
    REAL, DIMENSION(knon,klev+1), INTENT(IN) :: paprs 
    REAL, DIMENSION(knon,klev), INTENT(IN)   :: delp  
    REAL, DIMENSION(knon,klev), INTENT(IN)   :: temp
    REAL, DIMENSION(knon,klev), INTENT(IN)   :: qbs
    REAL, INTENT(IN)                         :: dtime

! Output arguments
!****************************************************************************************
    REAL, DIMENSION(knon), INTENT(OUT)       :: Acoef_QBS_out
    REAL, DIMENSION(knon), INTENT(OUT)       :: Bcoef_QBS_out

    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: Ccoef_QBS_out
    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: Dcoef_QBS_out
    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: Kcoef_qbs_out
    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: gama_qbs_out

! Local variables
!****************************************************************************************
    REAL :: yCcoef_QBS(knon,klev)
    REAL :: yDcoef_QBS(knon,klev)
    REAL :: yAcoef_QBS(knon), yBcoef_QBS(knon)
    REAL :: yKcoefqbs(knon,klev)
    REAL :: ygamaqbs(1:knon,2:klev)

    REAL, DIMENSION(knon)                    :: psref 
    REAL                                     :: delz, pkh
    INTEGER                                  :: k, i, j, ierr


!****************************************************************************************
! 2)
! Definition of the coeficient K 
!
!****************************************************************************************
    yKcoefqbs(:,:) = 0.0
    DO k = 2, klev
       DO i = 1, knon
          yKcoefqbs(i,k) = &
               coefqbs(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
       ygamaqbs(:,:) = 0.0
 
! 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 de contenu en neige soufflee en difference de neige soufflee entre centre de couches
             ygamaqbs(i,k) = ygamaqbs(i,k) * delz    
          ENDDO
       ENDDO

    ELSE
       ygamaqbs(:,:) = 0.0
    ENDIF
    

!****************************************************************************************    
! 4)
! Calculte the coefficients C and D for specific content of blowing snow, qbs
!
!****************************************************************************************
    
    CALL calc_coef_qbs(knon, yKcoefqbs(:,:), ygamaqbs(:,:), delp(:,:), qbs(:,:), &
         yCcoef_QBS(:,:), yDcoef_QBS(:,:), yAcoef_QBS, yBcoef_QBS)

 
!****************************************************************************************
! 5)
! Return the first layer in output variables
!
!****************************************************************************************
    Acoef_QBS_out = yAcoef_QBS
    Bcoef_QBS_out = yBcoef_QBS

!****************************************************************************************
! 6)
! If Pbl is split, return also the other layers in output variables
!
!****************************************************************************************
    IF (mod(iflag_pbl_split,10) .ge.1) THEN
    DO k= 1, klev
      DO i= 1, knon
        Ccoef_QBS_out(i,k) = yCcoef_QBS(i,k)
        Dcoef_QBS_out(i,k) = yDcoef_QBS(i,k)
        Kcoef_qbs_out(i,k) = yKcoefqbs(i,k)
          IF (k.eq.1) THEN
            gama_qbs_out(i,k)  = 0.
          ELSE
            gama_qbs_out(i,k)  = ygamaqbs(i,k)
          ENDIF
      ENDDO
    ENDDO
!!!      
       ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
!!!
    DO k= 1, klev
      DO j= 1, knon
        i=ni(j)
        IF (k==1) THEN
          Acoef_QBS(i) = yAcoef_QBS(j)
          Bcoef_QBS(i) = yBcoef_QBS(j)
        ENDIF
        IF (k>=2) gamaqbs(i,k)=ygamaqbs(j,k)
        Ccoef_QBS(i,k) = yCcoef_QBS(j,k)
        Dcoef_QBS(i,k) = yDcoef_QBS(j,k)
        Kcoefqbs(i,k) = yKcoefqbs(j,k)
      ENDDO
    ENDDO

  END SUBROUTINE climb_qbs_down
!
!****************************************************************************************
!
  SUBROUTINE calc_coef_qbs(knon, Kcoef, gama, delp, X, Ccoef, Dcoef, Acoef, Bcoef)
!$gpum horizontal knon
!
! Calculate the coefficients C and D in : X(k) = C(k) + D(k)*X(k-1)
! where X is QQBS, and k the vertical level k=1,klev
USE yomcst_mod_h
USE dimphy, ONLY : klev

! Input arguments
!****************************************************************************************
    INTEGER, INTENT(IN)                      :: knon
    REAL, DIMENSION(knon,klev), INTENT(IN)   :: Kcoef, delp
    REAL, DIMENSION(knon,klev), INTENT(IN)   :: X
    REAL, DIMENSION(knon,2:klev), INTENT(IN) :: gama

! Output arguments
!****************************************************************************************
    REAL, DIMENSION(knon), INTENT(OUT)       :: Acoef, Bcoef
    REAL, DIMENSION(knon,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_qbs
!
!****************************************************************************************
!
  SUBROUTINE climb_qbs_up(knon, ni, dtime, qbs_old, &
       flx_qbs1, paprs, pplay, &
       Acoef_QBS_in, Bcoef_QBS_in, &
       Ccoef_QBS_in, Dcoef_QBS_in, &
       Kcoef_qbs_in, gama_qbs_in, &
       flux_qbs, d_qbs)
!$gpum horizontal knon
! 
! This routine calculates the flux and tendency of the specific content of blowing snow qbs 
! The quantity qbs is calculated according to 
! X(k) = C(k) + D(k)*X(k-1) for X=[qbs], where the coefficients 
! C and D are known from before and k is index of the vertical layer.
!   
USE yomcst_mod_h
USE compbl_mod_h
USE dimphy, ONLY : klev

! Input arguments
!****************************************************************************************
    INTEGER, INTENT(IN)                      :: knon
    INTEGER, INTENT(IN)                      :: ni(knon)
    REAL, INTENT(IN)                         :: dtime
    REAL, DIMENSION(knon,klev), INTENT(IN)   :: qbs_old
    REAL, DIMENSION(knon), INTENT(IN)        :: flx_qbs1
    REAL, DIMENSION(knon,klev+1), INTENT(IN) :: paprs
    REAL, DIMENSION(knon,klev), INTENT(IN)   :: pplay

!!! nrlmd le 02/05/2011
    REAL, DIMENSION(knon), INTENT(IN)        :: Acoef_QBS_in, Bcoef_QBS_in
    REAL, DIMENSION(knon,klev), INTENT(IN)   :: Ccoef_QBS_in, Dcoef_QBS_in
    REAL, DIMENSION(knon,klev), INTENT(IN)   :: Kcoef_qbs_in, gama_qbs_in
!!!

! Output arguments
!****************************************************************************************
    REAL, DIMENSION(knon,klev), INTENT(OUT)  :: flux_qbs, d_qbs

! Local variables
!****************************************************************************************
    REAL :: yCcoef_QBS(knon,klev)
    REAL :: yDcoef_QBS(knon,klev)
    REAL :: yAcoef_QBS(knon), yBcoef_QBS(knon)
    REAL :: yKcoefqbs(knon,klev)
    REAL :: ygamaqbs(1:knon,2:klev)

    REAL, DIMENSION(knon,klev)               :: qbs_new
    REAL, DIMENSION(knon)                    :: psref         
    INTEGER                                  :: k, i, j, ierr
!****************************************************************************************
! 1) 
! Definition of some variables
    REAL, DIMENSION(knon,klev)               :: zairm
!
!****************************************************************************************
    flux_qbs(:,:) = 0.0
    d_qbs(:,:)    = 0.0

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

    DO k= 1, klev
      DO j= 1, knon
        i=ni(j)
        IF (k==1) THEN
          yAcoef_QBS(j) = Acoef_QBS(i)
          yBcoef_QBS(j) = Bcoef_QBS(i)
        ENDIF
        IF (k>=2) ygamaqbs(j,k)=gamaqbs(i,k)
        yCcoef_QBS(j,k) = Ccoef_QBS(i,k)
        yDcoef_QBS(j,k) = Dcoef_QBS(i,k)
        yKcoefqbs(j,k) = Kcoefqbs(i,k)
      ENDDO
    ENDDO

       IF (mod(iflag_pbl_split,10) .ge.1) THEN
    DO i = 1, knon
      yAcoef_QBS(i)=Acoef_QBS_in(i)
      yBcoef_QBS(i)=Bcoef_QBS_in(i)
    ENDDO
    DO k = 1, klev
      DO i = 1, knon
        yCcoef_QBS(i,k)=Ccoef_QBS_in(i,k)
        yDcoef_QBS(i,k)=Dcoef_QBS_in(i,k)
        yKcoefqbs(i,k)=Kcoef_qbs_in(i,k)
          IF (k.gt.1) THEN
            ygamaqbs(i,k)=gama_qbs_in(i,k)
          ENDIF
      ENDDO
    ENDDO
!!!      
       ENDIF  ! (mod(iflag_pbl_split,2) .ge.1)
!!!

!****************************************************************************************
! 2)
! Calculation of QBS
!
!****************************************************************************************

!- First layer
    qbs_new(1:knon,1) = yAcoef_QBS(1:knon) + yBcoef_QBS(1:knon)*flx_qbs1(1:knon)*dtime
!- All the other layers 
    DO k = 2, klev
       DO i = 1, knon
          qbs_new(i,k) = yCcoef_QBS(i,k) + yDcoef_QBS(i,k)*qbs_new(i,k-1)
       END DO
    END DO
!****************************************************************************************
! 3)
! Calculation of the flux for QBS
!
!****************************************************************************************

!- The flux at first layer, k=1
    flux_qbs(1:knon,1)=flx_qbs1(1:knon)

!- The flux at all layers above surface
    DO k = 2, klev
       DO i = 1, knon
          flux_qbs(i,k) = (yKcoefqbs(i,k)/RG/dtime) * &
               (qbs_new(i,k)-qbs_new(i,k-1)+ygamaqbs(i,k))
       END DO
    END DO

!****************************************************************************************
! 4)
! Calculation of tendency for QBS
!
!****************************************************************************************
    DO k = 1, klev
       DO i = 1, knon
          d_qbs(i,k) = qbs_new(i,k) - qbs_old(i,k)
          zairm(i, k) = (paprs(i,k)-paprs(i,k+1))/rg
        END DO
    END DO

    DO k= 1, klev
      DO j= 1, knon
        i=ni(j)
        IF (k==1) THEN
          Acoef_QBS(i) = yAcoef_QBS(j)
          Bcoef_QBS(i) = yBcoef_QBS(j)
        ENDIF
        IF (k>=2) gamaqbs(i,k)=ygamaqbs(j,k)
        Ccoef_QBS(i,k) = yCcoef_QBS(j,k)
        Dcoef_QBS(i,k) = yDcoef_QBS(j,k)
        Kcoefqbs(i,k) = yKcoefqbs(j,k)
      ENDDO
    ENDDO

  END SUBROUTINE climb_qbs_up
!
!****************************************************************************************
!
END MODULE climb_qbs_mod

 




