source: LMDZ6/trunk/libf/phylmd/lmdz_blosno_sublimsedim.f90

module lmdz_blosno_sublimsedim

contains

!==============================================================================
   subroutine blosno_sublimsedim(ngrid,nlay,dtime,ustar,temp,qv,qb,pplay,paprs,dtemp_bs,dqv_bs,dqb_bs,bsfl,precip_bs)

!==============================================================================
! Routine that calculates the sublimation, melting and sedimentation
! of blowing snow
! Reference: Vignon et al 2025, GMD, https://doi.org/10.5194/egusphere-2025-2871
!
! Etienne Vignon, October 2022
!==============================================================================

      use lmdz_blosno_ini, only: iflag_sublim_bs, iflag_sedim_bs, coef_sub_bs, RTT, RD, RG, fallv_bs
      use lmdz_blosno_ini, only: qbmin, RCPD, RLSTT, RLMLT, RLVTT, RVTMP2, RV, RPI
      use lmdz_blosno_ini, only: tbsmelt, taumeltbs0, rhobs, r_bs
      USE lmdz_lscp_tools, only: calc_qsat_ecmwf

      implicit none

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

!INPUT
!=====
      integer, intent(in)                     :: ngrid ! number of horizontal grid points
      integer, intent(in)                     :: nlay  ! number of vertical layers
      real, intent(in)                        :: dtime ! physics time step [s]
      real, intent(in), dimension(ngrid)      :: ustar ! surface friction velocity [m/s]
      real, intent(in), dimension(ngrid, nlay) :: temp  ! temperature of the air [K]
      real, intent(in), dimension(ngrid, nlay) :: qv    ! specific content of water [kg/kg]
      real, intent(in), dimension(ngrid, nlay) :: qb    ! specific content of blowing snow [kg/kg]
      real, intent(in), dimension(ngrid, nlay) :: pplay ! pressure at the middle of layers [Pa]
      real, intent(in), dimension(ngrid, nlay + 1) :: paprs ! pressure at the layer bottom interface [Pa]

! OUTPUT
!========
      real, intent(out), dimension(ngrid, nlay) :: dtemp_bs ! temperature tendency [K/s]
      real, intent(out), dimension(ngrid, nlay) :: dqv_bs   ! water vapor tendency [kg/kg/s]
      real, intent(out), dimension(ngrid, nlay) :: dqb_bs   ! blowing snow mass tendancy [kg/kg/s]
      real, intent(out), dimension(ngrid, nlay + 1) :: bsfl   ! vertical profile of blowing snow vertical flux [kg/m2/s]
      real, intent(out), dimension(ngrid)      :: precip_bs ! surface sedimentation flux of blowing snow [kg/s/m2]

! LOCAL
!======

      integer                                  :: k, i, n
      real                                     :: cpd, cpw, dqbsub, maxdqbsub, dqbmelt, zmair
      real                                     :: dqsedim, precbs, dqvmelt, zmelt, taumeltbs
      real                                     :: maxdqvmelt, rhoair, dz, dqbsedim
      real                                     :: delta_p, b_p, a_p, c_p, c_sub, qvsub
      real                                     :: p0, T0, Dv, Aprime, Bprime, Ka, radius0, radius
      real, dimension(ngrid)                   :: ztemp, zqv, zqb, zpres, qsi, dqsi, qsl, dqsl, qzero, sedim
      real, dimension(ngrid)                   :: zpaprsup, zpaprsdn, ztemp_up, zqb_up, zvelo
      real, dimension(ngrid, nlay)              :: temp_seri, qb_seri, qv_seri, zz

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

      qzero(:) = 0.
      dtemp_bs(:, :) = 0.
      dqv_bs(:, :) = 0.
      dqb_bs(:, :) = 0.
      zvelo(:) = 0.
      sedim(:) = 0.
      precip_bs(:) = 0.
      bsfl(:, :) = 0.
      zz(:, :) = 0.

      p0 = 101325.0    ! ref pressure

      DO k = 1, nlay
         DO i = 1, ngrid
            temp_seri(i, k) = temp(i, k)
            qv_seri(i, k) = qv(i, k)
            qb_seri(i, k) = qb(i, k)
            zz(i, k) = zz(i, k) + (paprs(i, k) - paprs(i, k + 1))/(pplay(i, k)/RD/temp(i, k)*RG)
         END DO
      END DO

! Sedimentation scheme
!----------------------

      IF (iflag_sedim_bs .GT. 0) THEN
! begin of top-down loop
         DO k = nlay, 1, -1

            DO i = 1, ngrid
               ztemp(i) = temp_seri(i, k)
               zqv(i) = qv_seri(i, k)
               zqb(i) = qb_seri(i, k)
               zpres(i) = pplay(i, k)
               zpaprsup(i) = paprs(i, k + 1)
               zpaprsdn(i) = paprs(i, k)
            END DO

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

               DO i = 1, ngrid
                  zmair = (zpaprsdn(i) - zpaprsup(i))/RG
                  ! RVTMP2=rcpv/rcpd-1
                  cpd = RCPD*(1.0 + RVTMP2*zqv(i))
                  cpw = RCPD*RVTMP2
                  ! zqb_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
                  zqb_up(i) = sedim(i)*dtime/zmair
                  ztemp_up(i) = temp_seri(i, k + 1)

                  ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
                  ztemp(i) = (ztemp_up(i)*zqb_up(i)*cpw + cpd*ztemp(i)) &
                             /(cpd + zqb_up(i)*cpw)
               END DO

            END IF

            DO i = 1, ngrid

               rhoair = zpres(i)/ztemp(i)/RD
               dz = (zpaprsdn(i) - zpaprsup(i))/(rhoair*RG)
               ! BS fall velocity assumed to be constant for now
               zvelo(i) = fallv_bs
               ! dqb/dt_sedim=1/rho(sedim/dz - rho w qb/dz)
               ! implicit resolution
               dqbsedim = (sedim(i)/rhoair/dz*dtime + zqb(i))/(1.+zvelo(i)*dtime/dz) - zqb(i)
               ! flux and dqb update
               zqb(i) = zqb(i) + dqbsedim
               !sedim(i) = sedim(i)-dqbsedim/dtime*rhoair*dz
               sedim(i) = rhoair*zvelo(i)*zqb(i)

               ! variables update:
               bsfl(i, k) = sedim(i)

               qb_seri(i, k) = zqb(i)
               qv_seri(i, k) = zqv(i)
               temp_seri(i, k) = ztemp(i)

            END DO ! Loop on ngrid

         END DO ! vertical loop

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

      END IF

!+++++++++++++++++++++++++++++++++++++++++++++++
! Sublimation and melting
!++++++++++++++++++++++++++++++++++++++++++++++
      IF (iflag_sublim_bs .GT. 0) THEN

         DO k = 1, nlay

            DO i = 1, ngrid
               ztemp(i) = temp_seri(i, k)
               zqv(i) = qv_seri(i, k)
               zqb(i) = qb_seri(i, k)
               zpres(i) = pplay(i, k)
               zpaprsup(i) = paprs(i, k + 1)
               zpaprsdn(i) = paprs(i, k)

               ! computation of blowing snow mean radius. Either constant value = r_bs if >0
               ! or use of the height-dependent formulation from Sauter et al. 2013 and Saigger et al. 2024
               IF (r_bs .GT. 0.) THEN
                  radius = r_bs
               ELSE
                  radius0 = 0.5*(ustar(i)*(7.8e-6)/0.036 + 31.e-6)
                  radius = radius0*zz(i, k)**(-0.258)
               END IF

            END DO

            ! calulation saturation specific humidity
            CALL CALC_QSAT_ECMWF(ngrid, ztemp(:), qzero(:), zpres(:), RTT, 2, .false., qsi(:), dqsi(:))

            DO i = 1, ngrid

               rhoair = zpres(i)/ztemp(i)/RD
               dz = (zpaprsdn(i) - zpaprsup(i))/(rhoair*RG)
               ! BS fall velocity assumed to be constant for now
               zvelo(i) = fallv_bs

               IF (ztemp(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., (ztemp(i) - RTT)/(tbsmelt - RTT)))
                  dqvmelt = min(zqb(i), -1.*zqb(i)*(exp(-dtime/taumeltbs) - 1.))
                  maxdqvmelt = max(RCPD*(1.0 + RVTMP2*(zqv(i) + zqb(i)))*(ztemp(i) - RTT)/(RLMLT + RLVTT), 0.)
                  dqvmelt = min(dqvmelt, maxdqvmelt)
                  ! qv update, melting + evaporation
                  zqv(i) = zqv(i) + dqvmelt
                  ! temp update melting + evaporation
                  ztemp(i) = ztemp(i) - dqvmelt*(RLMLT + RLVTT)/RCPD/(1.0 + RVTMP2*(zqv(i) + zqb(i)))
                  ! qb update melting + evaporation
                  zqb(i) = zqb(i) - dqvmelt

               ELSE
                  ! negative celcius temperature
                  ! Sublimation scheme

                  ! Sublimation formulation for ice crystals from Pruppacher & Klett, Rutledge & Hobbs 1983
                  ! assuming monodispered crystal distrib
                  ! dqb/dt_sub=-coef_sub_bs*(1-qv/qsi)*nc*8*radius/(Aprime+Bprime)
                  ! assuming Mi=rhobs*4/3*pi*radius**3
                  ! qb=nc*Mi -> nc=qb/Mi
                  ! dqb/dt_sub=-coef_sub_bs*(1-qv/qsi)*6*qb/(rhobs*pi*radius**2)/(Aprime+Bprime)
                  ! dqb/dt_sub=-c_sub(1-qv/qsi)*qb
                  ! c_sub=coef_sub_bs*6/(rhobs*pi*radius**2)/(Aprime+Bprime)
                  !
                  ! Note the strong coupling between specific contents of water vapor and blowing snow during sublimation
                  ! equations dqv/dt_sub and dqb/dt_sub must be solved jointly to prevent irrealistic increase of water vapor
                  ! at typical physics time steps
                  ! we thus solve the differential equation using an implicit scheme for both qb and qv

                  ! we do not consider deposition, only sublimation
                  IF (zqv(i) .LT. qsi(i)) THEN
                     rhoair = zpres(i)/ztemp(i)/RD
                     Dv = 0.0001*0.211*(p0/zpres(i))*((ztemp(i)/RTT)**1.94)           ! water vapor diffusivity in air, SI
                     Ka = (5.69 + 0.017*(ztemp(i) - RTT))*1.e-5*100.*4.184                ! thermal conductivity of the air, SI
                     Aprime = RLSTT/Ka/ztemp(i)*(RLSTT/RV/ztemp(i) - 1.)
                     Bprime = 1./(rhoair*Dv*qsi(i))
                     c_sub = coef_sub_bs*3./(rhobs*radius**2)/(Aprime + Bprime)
                     c_p = -zqb(i)
                     b_p = 1.+c_sub*dtime - c_sub*dtime/qsi(i)*zqb(i) - c_sub*dtime/qsi(i)*zqv(i)
                     a_p = c_sub*dtime/qsi(i)
                     delta_p = (b_p**2) - 4.*a_p*c_p
                     dqbsub = (-b_p + sqrt(delta_p))/(2.*a_p) - zqb(i)
                     dqbsub = MIN(0.0, MAX(dqbsub, -zqb(i)))
                     ! Sublimation limit: we ensure that the whole mesh does not exceed saturation wrt ice
                     maxdqbsub = MAX(0.0, qsi(i) - zqv(i))
                     dqbsub = MAX(dqbsub, -maxdqbsub)
                  ELSE
                     dqbsub = 0.
                  END IF

                  ! vapor, temperature, precip fluxes update following sublimation
                  zqv(i) = zqv(i) - dqbsub
                  zqb(i) = zqb(i) + dqbsub
                  ztemp(i) = ztemp(i) + dqbsub*RLSTT/RCPD/(1.0 + RVTMP2*(zqv(i) + zqb(i)))

               END IF

               ! if qb<qbmin, sublimate or melt and evaporate qb
               ! see Gerber et al. 2023, JGR Atmos for the choice of qbmin

               IF (zqb(i) .LT. qbmin) THEN
                  zqv(i) = zqv(i) + zqb(i)
                  IF (ztemp(i) .LT. RTT) THEN
                     ztemp(i) = ztemp(i) - zqb(i)*RLSTT/RCPD/(1.0 + RVTMP2*(zqv(i)))
                  ELSE
                     ztemp(i) = ztemp(i) - zqb(i)*(RLVTT + RLMLT)/RCPD/(1.0 + RVTMP2*(zqv(i)))
                  END IF
                  zqb(i) = 0.
               END IF

               ! variables update
               temp_seri(i, k) = ztemp(i)
               qv_seri(i, k) = zqv(i)
               qb_seri(i, k) = zqb(i)
            END DO
         END DO

      END IF

! OUTPUTS
!++++++++++

! 3D variables
      DO k = 1, nlay
         DO i = 1, ngrid
            dqb_bs(i, k) = qb_seri(i, k) - qb(i, k)
            dqv_bs(i, k) = qv_seri(i, k) - qv(i, k)
            dtemp_bs(i, k) = temp_seri(i, k) - temp(i, k)
         END DO
      END DO

      return

   end subroutine blosno_sublimsedim
!==============================================================================

end module lmdz_blosno_sublimsedim