source: LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_sublim_sedim.F90 @ 4821

module lmdz_blowing_snow_sublim_sedim

contains
subroutine blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,qv,qbs,pplay,paprs,dtemp_bs,dq_bs,dqbs_bs,bsfl,precip_bs)

!==============================================================================
! Routine that calculates the evaporation and sedimentation of blowing snow
! inspired by what is done in lscp_mod
! Etienne Vignon, October 2022
!==============================================================================


use lmdz_blowing_snow_ini, only : coef_eva_bs,RTT,RD,RG,expo_eva_bs, fallv_bs, qbsmin
use lmdz_blowing_snow_ini, only : RCPD, RLSTT, RLMLT, RLVTT, RVTMP2, tbsmelt, taumeltbs0
USE lmdz_lscp_tools, only : calc_qsat_ecmwf

implicit none


!++++++++++++++++++++++++++++++++++++++++++++++++++++
! Declarations
!++++++++++++++++++++++++++++++++++++++++++++++++++++

!INPUT
!=====

integer, intent(in)                     :: ngrid,nlay
real, intent(in)                        :: dtime
real, intent(in), dimension(ngrid,nlay) :: temp
real, intent(in), dimension(ngrid,nlay) :: qv
real, intent(in), dimension(ngrid,nlay) :: qbs
real, intent(in), dimension(ngrid,nlay) :: pplay
real, intent(in), dimension(ngrid,nlay+1) :: paprs



! OUTPUT
!========


real, intent(out), dimension(ngrid,nlay) :: dtemp_bs
real, intent(out), dimension(ngrid,nlay) :: dq_bs
real, intent(out), dimension(ngrid,nlay) :: dqbs_bs
real, intent(out), dimension(ngrid,nlay+1) :: bsfl
real, intent(out), dimension(ngrid)      :: precip_bs


! LOCAL
!======


integer                                  :: k,i,n
real                                     :: cpd, cpw, dqsub, maxdqsub, dqbsmelt
real                                     :: dqsedim,precbs, dqmelt, zmelt, taumeltbs
real                                     :: maxdqmelt, rhoair, dz
real, dimension(ngrid)                   :: zt,zq,zqbs,qsi,dqsi,qsl, dqsl,qzero,sedim,sedimn
real, dimension(ngrid)                   :: zqbsi, zqbs_up, zmair
real, dimension(ngrid)                   :: zvelo 

!++++++++++++++++++++++++++++++++++++++++++++++++++
! Initialisation
!++++++++++++++++++++++++++++++++++++++++++++++++++

qzero(:)=0.
dtemp_bs(:,:)=0.
dq_bs(:,:)=0.
dqbs_bs(:,:)=0.
zvelo(:)=0.
zt(:)=0.
zq(:)=0.
zqbs(:)=0.
sedim(:)=0.
sedimn(:)=0.

! begin of top-down loop
DO k = nlay, 1, -1
    
    DO i=1,ngrid
        zt(i)=temp(i,k)
        zq(i)=qv(i,k)
        zqbs(i)=qbs(i,k)
    ENDDO

     ! thermalization of blowing snow precip coming from above     
    IF (k.LE.nlay-1) THEN

        DO i = 1, ngrid
            zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
            ! RVTMP2=rcpv/rcpd-1
            cpd=RCPD*(1.0+RVTMP2*zq(i))
            cpw=RCPD*RVTMP2
            ! zqbs_up: blowing snow mass that has to be thermalized with 
            ! layer's air so that precipitation at the ground has the
            ! same temperature as the lowermost layer
            zqbs_up(i) = sedim(i)*dtime/zmair(i)
            ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
            zt(i) = ( (temp(i,k+1)+dtemp_bs(i,k+1))*zqbs_up(i)*cpw + cpd*zt(i) ) &
             / (cpd + zqbs_up(i)*cpw)
        ENDDO
    ENDIF



    ! calulation saturation specific humidity
    CALL CALC_QSAT_ECMWF(ngrid,zt(:),qzero(:),pplay(:,k),RTT,2,.false.,qsi(:),dqsi(:))
    CALL CALC_QSAT_ECMWF(ngrid,zt(:),qzero(:),pplay(:,k),RTT,1,.false.,qsl(:),dqsl(:))
    
    
    ! 3 processes: melting, sublimation and precipitation of blowing snow
    DO i = 1, ngrid

          IF (zt(i) .GT. RTT) THEN
             
             ! if temperature is positive, we assume that part of the blowing snow 
             ! already present  melts and evaporates with a typical time 
             ! constant taumeltbs
           
             taumeltbs=taumeltbs0*exp(-max(0.,(zt(i)-RTT)/(tbsmelt-RTT)))
             dqmelt=min(zqbs(i),zqbs(i)/taumeltbs*dtime)
             maxdqmelt= max(RCPD*(1.0+RVTMP2*(zq(i)+zqbs(i)))*(zt(i)-RTT)/(RLMLT+RLVTT),0.)
             dqmelt=min(max(dqmelt,0.),maxdqmelt)
             ! q update, melting + evaporation
             zq(i) = zq(i) + dqmelt
             ! temp update melting + evaporation
             zt(i) = zt(i) - dqmelt * (RLMLT+RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
             ! qbs update melting + evaporation
             zqbs(i)=zqbs(i)-dqmelt

          ELSE 
              ! negative celcius temperature     
              ! Sublimation scheme   

              rhoair=pplay(i,k)/zt(i)/RD 
              dqsub = 1./rhoair*coef_eva_bs*(1.0-zq(i)/qsi(i))*((rhoair*zvelo(i)*zqbs(i))**expo_eva_bs)*dtime
              dqsub = MAX(0.0,MIN(dqsub,zqbs(i)))

              ! Sublimation limit: we ensure that the whole mesh does not exceed saturation wrt ice
              maxdqsub = MAX(0.0, qsi(i)-zq(i))
              dqsub = MIN(dqsub,maxdqsub)


              ! vapor, temperature, precip fluxes update following sublimation
              zq(i) = zq(i) + dqsub     
              zqbs(i)=zqbs(i)-dqsub           
              zt(i) = zt(i) - dqsub*RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))

          ENDIF

        ! Sedimentation scheme
           
        rhoair  = pplay(i,k) / zt(i) / RD
        dz      = (paprs(i,k)-paprs(i,k+1)) / (rhoair*RG)
        ! BS fall velocity assumed to be constant for now
        zvelo(i) = fallv_bs

        sedimn(i) = rhoair*zqbs(i)*zvelo(i)

        ! exact numerical resolution of sedimentation
        ! assuming fall velocity is constant 

        zqbs(i) = zqbs(i)*exp(-zvelo(i)/dz*dtime)+sedim(i)/rhoair/zvelo(i)

        ! flux update
        sedim(i) = sedimn(i)



        ! if qbs<qbsmin, sublimate or melt and evaporate qbs
        ! see Gerber et al. 2023, JGR Atmos for the choice of qbsmin
        IF (zqbs(i) .LT. qbsmin) THEN                       
              zq(i) = zq(i)+zqbs(i)
              IF (zt(i) .LT. RTT) THEN
                 zt(i) = zt(i) - zqbs(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
              ELSE
               zt(i) = zt(i) - zqbs(i) * (RLVTT+RLMLT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
              ENDIF
              zqbs(i)=0.
        ENDIF


    ! Outputs:
        bsfl(i,k)=sedim(i)
        dqbs_bs(i,k) = zqbs(i)-qbs(i,k)
        dq_bs(i,k) = zq(i) - qv(i,k)
        dtemp_bs(i,k) = zt(i) - temp(i,k)

    ENDDO ! Loop on ngrid


ENDDO ! vertical loop


!surface bs flux
DO i = 1, ngrid
  precip_bs(i) = sedim(i)
ENDDO


return

end subroutine blowing_snow_sublim_sedim
end module lmdz_blowing_snow_sublim_sedim