module lmdz_blowing_snow_sublim_sedim contains subroutine blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,qv,qb,pplay,paprs,dtemp_bs,dqv_bs,dqb_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 : iflag_sublim_bs, iflag_sedim_bs, coef_sub_bs,RTT,RD,RG,fallv_bs use lmdz_blowing_snow_ini, only : qbmin, RCPD, RLSTT, RLMLT, RLVTT, RVTMP2, RV, RPI, tbsmelt, taumeltbs0, rhobs, r_bs 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) :: qb 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) :: dqv_bs real, intent(out), dimension(ngrid,nlay) :: dqb_bs real, intent(out), dimension(ngrid,nlay+1) :: bsfl real, intent(out), dimension(ngrid) :: precip_bs ! 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 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 !++++++++++++++++++++++++++++++++++++++++++++++++++ ! Initialisation !++++++++++++++++++++++++++++++++++++++++++++++++++ qzero(:)=0. dtemp_bs(:,:)=0. dqv_bs(:,:)=0. dqb_bs(:,:)=0. zvelo(:)=0. sedim(:)=0. precip_bs(:)=0. bsfl(:,:)=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) ENDDO ENDDO ! 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) ENDDO ! 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) ENDDO ENDIF 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) ENDDO ! Loop on ngrid ENDDO ! vertical loop !surface bs flux DO i = 1, ngrid precip_bs(i) = sedim(i) ENDDO ENDIF !+++++++++++++++++++++++++++++++++++++++++++++++ ! 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) ENDDO ! 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*r_bs/(Aprime+Bprime) ! assuming Mi=rhobs*4/3*pi*r_bs**3 ! rhoair qb=nc*Mi -> nc=rhoair qb/Mi ! dqb/dt_sub=-coef_sub_bs*(1-qv/qsi)*6*rhoair*qb/(rhobs*pi*r_bs**2)/(Aprime+Bprime) ! dqb/dt_sub=-c_sub(1-qv/qsi)*qb ! c_sub=coef_sub_bs*6*rhoair/(rhobs*pi*r_bs**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*6.*rhoair/(rhobs*RPI*r_bs**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. ENDIF ! 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))) ENDIF ! 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))) ENDIF zqb(i)=0. ENDIF ! variables update temp_seri(i,k)=ztemp(i) qv_seri(i,k)=zqv(i) qb_seri(i,k)=zqb(i) ENDDO ENDDO ENDIF ! 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) ENDDO ENDDO return end subroutine blowing_snow_sublim_sedim end module lmdz_blowing_snow_sublim_sedim