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 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 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),-1.*zqbs(i)*(exp(-dtime/taumeltbs)-1.)) maxdqmelt= max(RCPD*(1.0+RVTMP2*(zq(i)+zqbs(i)))*(zt(i)-RTT)/(RLMLT+RLVTT),0.) dqmelt=min(dqmelt,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 !---------------------- 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