1 |
|
|
subroutine blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,q,qbs,pplay,paprs,dtemp_bs,dq_bs,dqbs_bs,bsfl,precip_bs) |
2 |
|
|
|
3 |
|
|
!============================================================================== |
4 |
|
|
! Routine that calculates the evaporation and sedimentation of blowing snow |
5 |
|
|
! inspired by what is done in lscp_mod |
6 |
|
|
! Etienne Vignon, October 2022 |
7 |
|
|
!============================================================================== |
8 |
|
|
|
9 |
|
|
|
10 |
|
|
use blowing_snow_ini_mod, only : coef_eva_bs,RTT,RD,RG,expo_eva_bs, fallv_bs |
11 |
|
|
use blowing_snow_ini_mod, only : RCPD, RLSTT, RLMLT, RLVTT, RVTMP2, niter_bs |
12 |
|
|
USE lscp_tools_mod, only : calc_qsat_ecmwf |
13 |
|
|
|
14 |
|
|
implicit none |
15 |
|
|
|
16 |
|
|
|
17 |
|
|
!++++++++++++++++++++++++++++++++++++++++++++++++++++ |
18 |
|
|
! Declarations |
19 |
|
|
!++++++++++++++++++++++++++++++++++++++++++++++++++++ |
20 |
|
|
|
21 |
|
|
!INPUT |
22 |
|
|
!===== |
23 |
|
|
|
24 |
|
|
integer, intent(in) :: ngrid,nlay |
25 |
|
|
real, intent(in) :: dtime |
26 |
|
|
real, intent(in), dimension(ngrid,nlay) :: temp |
27 |
|
|
real, intent(in), dimension(ngrid,nlay) :: q |
28 |
|
|
real, intent(in), dimension(ngrid,nlay) :: qbs |
29 |
|
|
real, intent(in), dimension(ngrid,nlay) :: pplay |
30 |
|
|
real, intent(in), dimension(ngrid,nlay+1) :: paprs |
31 |
|
|
|
32 |
|
|
|
33 |
|
|
|
34 |
|
|
! OUTPUT |
35 |
|
|
!======== |
36 |
|
|
|
37 |
|
|
|
38 |
|
|
real, intent(out), dimension(ngrid,nlay) :: dtemp_bs |
39 |
|
|
real, intent(out), dimension(ngrid,nlay) :: dq_bs |
40 |
|
|
real, intent(out), dimension(ngrid,nlay) :: dqbs_bs |
41 |
|
|
real, intent(out), dimension(ngrid,nlay+1) :: bsfl |
42 |
|
|
real, intent(out), dimension(ngrid) :: precip_bs |
43 |
|
|
|
44 |
|
|
|
45 |
|
|
! LOCAL |
46 |
|
|
!====== |
47 |
|
|
|
48 |
|
|
|
49 |
|
|
integer :: k,i,n |
50 |
|
|
real :: zqev0, zqevi, zqevti, zcpair, zcpeau, dqbsmelt |
51 |
|
|
real :: dqsedim,precbs |
52 |
|
|
real, dimension(ngrid) :: zt,zq,zqbs,qsi,dqsi,qsl, dqsl,qzero,sedim,sedimn |
53 |
|
|
real, dimension(ngrid) :: zqbsi,zmqc, zmair, zdz |
54 |
|
|
real, dimension(ngrid,nlay) :: velo, zrho |
55 |
|
|
|
56 |
|
|
!++++++++++++++++++++++++++++++++++++++++++++++++++ |
57 |
|
|
! Initialisation |
58 |
|
|
!++++++++++++++++++++++++++++++++++++++++++++++++++ |
59 |
|
|
|
60 |
|
|
qzero(:)=0. |
61 |
|
|
dtemp_bs(:,:)=0. |
62 |
|
|
dq_bs(:,:)=0. |
63 |
|
|
dqbs_bs(:,:)=0. |
64 |
|
|
velo(:,:)=0. |
65 |
|
|
zt(:)=0. |
66 |
|
|
zq(:)=0. |
67 |
|
|
zqbs(:)=0. |
68 |
|
|
sedim(:)=0. |
69 |
|
|
|
70 |
|
|
! begin of top-down loop |
71 |
|
|
DO k = nlay, 1, -1 |
72 |
|
|
|
73 |
|
|
DO i=1,ngrid |
74 |
|
|
zt(i)=temp(i,k) |
75 |
|
|
zq(i)=q(i,k) |
76 |
|
|
zqbs(i)=qbs(i,k) |
77 |
|
|
ENDDO |
78 |
|
|
|
79 |
|
|
|
80 |
|
|
IF (k.LE.nlay-1) THEN |
81 |
|
|
|
82 |
|
|
! thermalization of blowing snow precip coming from above |
83 |
|
|
DO i = 1, ngrid |
84 |
|
|
zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG |
85 |
|
|
! RVTMP2=rcpv/rcpd-1 |
86 |
|
|
zcpair=RCPD*(1.0+RVTMP2*zq(i)) |
87 |
|
|
zcpeau=RCPD*RVTMP2 |
88 |
|
|
! zmqc: precipitation mass that has to be thermalized with |
89 |
|
|
! layer's air so that precipitation at the ground has the |
90 |
|
|
! same temperature as the lowermost layer |
91 |
|
|
zmqc(i) = (sedim(i))*dtime/zmair(i) |
92 |
|
|
! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer |
93 |
|
|
zt(i) = ( (temp(i,k+1)+dtemp_bs(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) & |
94 |
|
|
/ (zcpair + zmqc(i)*zcpeau) |
95 |
|
|
ENDDO |
96 |
|
|
ELSE |
97 |
|
|
DO i = 1, ngrid |
98 |
|
|
zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG |
99 |
|
|
zmqc(i) = 0. |
100 |
|
|
ENDDO |
101 |
|
|
|
102 |
|
|
ENDIF |
103 |
|
|
|
104 |
|
|
|
105 |
|
|
|
106 |
|
|
! calulation saturation specific humidity |
107 |
|
|
CALL CALC_QSAT_ECMWF(ngrid,zt(:),qzero(:),pplay(:,k),RTT,2,.false.,qsi(:),dqsi(:)) |
108 |
|
|
CALL CALC_QSAT_ECMWF(ngrid,zt(:),qzero(:),pplay(:,k),RTT,1,.false.,qsl(:),dqsl(:)) |
109 |
|
|
! sublimation calculation |
110 |
|
|
! SUndqvist formula dP/dz=beta*(1-q/qsat)*sqrt(P) |
111 |
|
|
|
112 |
|
|
DO i = 1, ngrid |
113 |
|
|
! if sedimentation: |
114 |
|
|
IF (sedim(i) .GT. 0.) THEN |
115 |
|
|
|
116 |
|
|
IF (zt(i) .GT. RTT) THEN |
117 |
|
|
! if positive celcius temperature, we assume |
118 |
|
|
! that all the blowing snow melt and evaporate |
119 |
|
|
zqevti=sedim(i)*(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG |
120 |
|
|
! we ensure that the whole mesh does not exceed saturation wrt liquid |
121 |
|
|
zqev0 = MAX(0.0, qsl(i)-zq(i)) |
122 |
|
|
zqevi = MIN(zqev0,zqevti) |
123 |
|
|
! New solid precipitation fluxes |
124 |
|
|
sedimn(i) = Max(0.,sedim(i) - zqevi*(paprs(i,k)-paprs(i,k+1))/RG/dtime) |
125 |
|
|
! vapor, temperature, precip fluxes update |
126 |
|
|
zq(i) = zq(i) - (sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime |
127 |
|
|
zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime |
128 |
|
|
! melting |
129 |
|
|
zt(i) = zt(i) & |
130 |
|
|
+ (sedimn(i)-sedim(i)) & |
131 |
|
|
* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime & |
132 |
|
|
* RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i))) |
133 |
|
|
! evaporation |
134 |
|
|
zt(i) = zt(i) & |
135 |
|
|
+ (sedimn(i)-sedim(i)) & |
136 |
|
|
* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime & |
137 |
|
|
* RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i))) |
138 |
|
|
|
139 |
|
|
|
140 |
|
|
ELSE |
141 |
|
|
zqevti = coef_eva_bs*(1.0-zq(i)/qsi(i))*(sedim(i)**expo_eva_bs) & |
142 |
|
|
*(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG |
143 |
|
|
zqevti = MAX(0.0,MIN(zqevti,sedim(i)))*RG*dtime/(paprs(i,k)-paprs(i,k+1)) |
144 |
|
|
|
145 |
|
|
! Sublimation limit: we ensure that the whole mesh does not exceed saturation wrt ice |
146 |
|
|
zqev0 = MAX(0.0, qsi(i)-zq(i)) |
147 |
|
|
zqevi = MIN(zqev0,zqevti) |
148 |
|
|
|
149 |
|
|
! New solid precipitation fluxes |
150 |
|
|
sedimn(i) = Max(0.,sedim(i) - zqevi*(paprs(i,k)-paprs(i,k+1))/RG/dtime) |
151 |
|
|
|
152 |
|
|
|
153 |
|
|
! vapor, temperature, precip fluxes update |
154 |
|
|
zq(i) = zq(i) - (sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime |
155 |
|
|
zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime |
156 |
|
|
zt(i) = zt(i) & |
157 |
|
|
+ (sedimn(i)-sedim(i)) & |
158 |
|
|
* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime & |
159 |
|
|
* RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i))) |
160 |
|
|
ENDIF |
161 |
|
|
|
162 |
|
|
|
163 |
|
|
! sedim update |
164 |
|
|
sedim(i)=sedimn(i) |
165 |
|
|
|
166 |
|
|
|
167 |
|
|
ELSE |
168 |
|
|
sedim(i)=0. |
169 |
|
|
ENDIF ! if sedim > 0 |
170 |
|
|
|
171 |
|
|
zqbsi(i)=zqbs(i) |
172 |
|
|
|
173 |
|
|
ENDDO ! loop on ngrid |
174 |
|
|
|
175 |
|
|
! Now sedimention scheme (alike that in lscp) |
176 |
|
|
DO n = 1, niter_bs |
177 |
|
|
DO i = 1, ngrid |
178 |
|
|
zrho(i,k) = pplay(i,k) / zt(i) / RD |
179 |
|
|
zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG) |
180 |
|
|
velo(i,k) = fallv_bs |
181 |
|
|
dqsedim = dtime/REAL(niter_bs)/zdz(i)*zqbs(i)*velo(i,k) ! dqice/dt=1/rho*d(rho*wice*qice)/dz |
182 |
|
|
precbs = MIN(MAX(dqsedim,0.0),zqbs(i)) |
183 |
|
|
zqbs(i) = MAX(zqbs(i)-1.*precbs , 0.0) |
184 |
|
|
ENDDO !loop on ngrid |
185 |
|
|
ENDDO ! loop on niter_bs |
186 |
|
|
|
187 |
|
|
! add to non sublimated precip |
188 |
|
|
DO i=1,ngrid |
189 |
|
|
sedim(i) = sedim(i)+max(0.,zqbsi(i)-zqbs(i))*(paprs(i,k)-paprs(i,k+1))/(RG*dtime) |
190 |
|
|
ENDDO |
191 |
|
|
|
192 |
|
|
! Outputs: |
193 |
|
|
DO i = 1, ngrid |
194 |
|
|
bsfl(i,k)=sedim(i) |
195 |
|
|
dqbs_bs(i,k) = zqbs(i)-qbs(i,k) |
196 |
|
|
dq_bs(i,k) = zq(i) - q(i,k) |
197 |
|
|
dtemp_bs(i,k) = zt(i) - temp(i,k) |
198 |
|
|
ENDDO |
199 |
|
|
|
200 |
|
|
|
201 |
|
|
ENDDO ! vertical loop |
202 |
|
|
|
203 |
|
|
|
204 |
|
|
!surface bs flux |
205 |
|
|
DO i = 1, ngrid |
206 |
|
|
precip_bs(i) = sedim(i) |
207 |
|
|
ENDDO |
208 |
|
|
|
209 |
|
|
|
210 |
|
|
return |
211 |
|
|
|
212 |
|
|
end subroutine blowing_snow_sublim_sedim |