Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

lmdz_blowing_snow_sublim_sedim.F90 @ 4727

Last change on this file since 4727 was 4727, checked in by idelkadi, 8 months ago
Merged trunk changes -r4488:4726 LMDZ_ECRad branch
File size: 8.3 KB

Line
1	module lmdz_blowing_snow_sublim_sedim
2
3	contains
4	subroutine blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,qv,qbs,pplay,paprs,dtemp_bs,dq_bs,dqbs_bs,bsfl,precip_bs)
5
6	!==============================================================================
7	! Routine that calculates the evaporation and sedimentation of blowing snow
8	! inspired by what is done in lscp_mod
9	! Etienne Vignon, October 2022
10	!==============================================================================
11
12
13	use lmdz_blowing_snow_ini, only : coef_eva_bs,RTT,RD,RG,expo_eva_bs, fallv_bs, qbsmin
14	use lmdz_blowing_snow_ini, only : RCPD, RLSTT, RLMLT, RLVTT, RVTMP2, tbsmelt, taumeltbs0
15	USE lmdz_lscp_tools, only : calc_qsat_ecmwf
16
17	implicit none
18
19
20	!++++++++++++++++++++++++++++++++++++++++++++++++++++
21	! Declarations
22	!++++++++++++++++++++++++++++++++++++++++++++++++++++
23
24	!INPUT
25	!=====
26
27	integer, intent(in) :: ngrid,nlay
28	real, intent(in) :: dtime
29	real, intent(in), dimension(ngrid,nlay) :: temp
30	real, intent(in), dimension(ngrid,nlay) :: qv
31	real, intent(in), dimension(ngrid,nlay) :: qbs
32	real, intent(in), dimension(ngrid,nlay) :: pplay
33	real, intent(in), dimension(ngrid,nlay+1) :: paprs
34
35
36
37	! OUTPUT
38	!========
39
40
41	real, intent(out), dimension(ngrid,nlay) :: dtemp_bs
42	real, intent(out), dimension(ngrid,nlay) :: dq_bs
43	real, intent(out), dimension(ngrid,nlay) :: dqbs_bs
44	real, intent(out), dimension(ngrid,nlay+1) :: bsfl
45	real, intent(out), dimension(ngrid) :: precip_bs
46
47
48	! LOCAL
49	!======
50
51
52	integer :: k,i,n
53	real :: zqev0, zqevi, zqevti, zcpair, zcpeau, dqbsmelt
54	real :: dqsedim,precbs, deltaqchaud, zmelt, taumeltbs
55	real :: maxdeltaqchaud
56	real, dimension(ngrid) :: zt,zq,zqbs,qsi,dqsi,qsl, dqsl,qzero,sedim,sedimn
57	real, dimension(ngrid) :: zqbsi,zmqc, zmair, zdz
58	real, dimension(ngrid,nlay) :: velo, zrho
59
60	!++++++++++++++++++++++++++++++++++++++++++++++++++
61	! Initialisation
62	!++++++++++++++++++++++++++++++++++++++++++++++++++
63
64	qzero(:)=0.
65	dtemp_bs(:,:)=0.
66	dq_bs(:,:)=0.
67	dqbs_bs(:,:)=0.
68	velo(:,:)=0.
69	zt(:)=0.
70	zq(:)=0.
71	zqbs(:)=0.
72	sedim(:)=0.
73
74	! begin of top-down loop
75	DO k = nlay, 1, -1
76
77	DO i=1,ngrid
78	zt(i)=temp(i,k)
79	zq(i)=qv(i,k)
80	zqbs(i)=qbs(i,k)
81	ENDDO
82
83
84	IF (k.LE.nlay-1) THEN
85
86	! thermalization of blowing snow precip coming from above
87	DO i = 1, ngrid
88	zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
89	! RVTMP2=rcpv/rcpd-1
90	zcpair=RCPD(1.0+RVTMP2zq(i))
91	zcpeau=RCPD*RVTMP2
92	! zmqc: precipitation mass that has to be thermalized with
93	! layer's air so that precipitation at the ground has the
94	! same temperature as the lowermost layer
95	zmqc(i) = (sedim(i))*dtime/zmair(i)
96	! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
97	zt(i) = ( (temp(i,k+1)+dtemp_bs(i,k+1))zmqc(i)zcpeau + zcpair*zt(i) ) &
98	/ (zcpair + zmqc(i)*zcpeau)
99	ENDDO
100	ELSE
101	DO i = 1, ngrid
102	zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
103	zmqc(i) = 0.
104	ENDDO
105
106	ENDIF
107
108
109
110	! calulation saturation specific humidity
111	CALL CALC_QSAT_ECMWF(ngrid,zt(:),qzero(:),pplay(:,k),RTT,2,.false.,qsi(:),dqsi(:))
112	CALL CALC_QSAT_ECMWF(ngrid,zt(:),qzero(:),pplay(:,k),RTT,1,.false.,qsl(:),dqsl(:))
113	! sublimation calculation
114	! SUndqvist formula dP/dz=beta(1-q/qsat)sqrt(P)
115
116	DO i = 1, ngrid
117
118	zrho(i,k) = pplay(i,k) / zt(i) / RD
119	zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
120	! BS fall velocity
121	velo(i,k) = fallv_bs
122
123
124	IF (zt(i) .GT. RTT) THEN
125
126	! if positive celcius temperature, we assume
127	! that part of the the blowing snow flux melts and evaporates
128
129	! vapor, bs, temperature, precip fluxes update
130	zmelt = ((zt(i)-RTT)/(tbsmelt-RTT))
131	zmelt = MIN(MAX(zmelt,0.),1.)
132	sedimn(i)=sedim(i)*zmelt
133	deltaqchaud=-(sedimn(i)-sedim(i))(RG/(paprs(i,k)-paprs(i,k+1)))dtime
134	! max evap such as celcius temperature cannot become negative
135	maxdeltaqchaud= max(RCPD(1.0+RVTMP2(zq(i)+zqbs(i)))*(zt(i)-RTT)/(RLMLT+RLVTT),0.)
136
137	deltaqchaud=min(max(deltaqchaud,0.),maxdeltaqchaud)
138	zq(i) = zq(i) + deltaqchaud
139
140	! melting + evaporation
141	zt(i) = zt(i) - deltaqchaud * (RLMLT+RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
142
143	sedim(i)=sedimn(i)
144
145	! if temperature still positive, we assume that part of the blowing snow
146	! already present in the mesh melts and evaporates with a typical time
147	! constant between taumeltbs0 and 0
148
149	IF ( zt(i) .GT. RTT ) THEN
150	taumeltbs=taumeltbs0*exp(-max(0.,(zt(i)-RTT)/(tbsmelt-RTT)))
151	deltaqchaud=min(zqbs(i),zqbs(i)/taumeltbs*dtime)
152	maxdeltaqchaud= max(RCPD(1.0+RVTMP2(zq(i)+zqbs(i)))*(zt(i)-RTT)/(RLMLT+RLVTT),0.)
153	deltaqchaud=min(max(deltaqchaud,0.),maxdeltaqchaud)
154	zq(i) = zq(i) + deltaqchaud
155
156	! melting + evaporation
157	zt(i) = zt(i) - deltaqchaud * (RLMLT+RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
158	! qbs update
159	zqbs(i)=max(zqbs(i)-deltaqchaud,0.)
160	ENDIF
161
162	! now sedimentation scheme with an exact numerical resolution
163	! (correct if fall velocity is constant)
164	zqbsi(i)=zqbs(i)
165	zqbs(i) = zqbsi(i)exp(-velo(i,k)/zdz(i)dtime)+sedim(i)/zrho(i,k)/velo(i,k)
166	zqbs(i) = max(zqbs(i),0.)
167
168	! flux update
169	sedim(i) = sedim(i) + zrho(i,k)zdz(i)/dtime(zqbsi(i)-zqbs(i))
170	sedim(i) = max(0.,sedim(i))
171
172	ELSE
173	! negative celcius temperature
174	! Sublimation scheme
175
176	zqevti = coef_eva_bs(1.0-zq(i)/qsi(i))(sedim(i)**expo_eva_bs) &
177	(paprs(i,k)-paprs(i,k+1))/pplay(i,k)zt(i)*RD/RG
178	zqevti = MAX(0.0,MIN(zqevti,sedim(i)))RGdtime/(paprs(i,k)-paprs(i,k+1))
179
180	! Sublimation limit: we ensure that the whole mesh does not exceed saturation wrt ice
181	zqev0 = MAX(0.0, qsi(i)-zq(i))
182	zqevi = MIN(zqev0,zqevti)
183
184	! New solid precipitation fluxes
185	sedimn(i) = Max(0.,sedim(i) - zqevi*(paprs(i,k)-paprs(i,k+1))/RG/dtime)
186
187
188	! vapor, temperature, precip fluxes update following sublimation
189	zq(i) = zq(i) - (sedimn(i)-sedim(i))(RG/(paprs(i,k)-paprs(i,k+1)))dtime
190	zq(i) = max(0., zq(i))
191	zt(i) = zt(i) &
192	+ (sedimn(i)-sedim(i)) &
193	* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime &
194	* RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
195
196	sedim(i)=sedimn(i)
197	zqbsi(i)=zqbs(i)
198
199	! now sedimentation scheme with an exact numerical resolution
200	! (correct if fall velocity is constant)
201
202	zqbs(i) = zqbsi(i)exp(-velo(i,k)/zdz(i)dtime)+sedim(i)/zrho(i,k)/velo(i,k)
203	zqbs(i) = max(zqbs(i),0.)
204
205	! flux update
206	sedim(i) = sedim(i) + zrho(i,k)zdz(i)/dtime(zqbsi(i)-zqbs(i))
207	sedim(i) = max(0.,sedim(i))
208
209
210	ENDIF
211
212
213	! if qbs<qbsmin, sublimate or melt and evaporate qbs
214	! see Gerber et al. 2023, JGR Atmos for the choice of qbsmin
215	IF (zqbs(i) .LT. qbsmin) THEN
216	zq(i) = zq(i)+zqbs(i)
217	IF (zt(i) .LT. RTT) THEN
218	zt(i) = zt(i) - zqbs(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
219	ELSE
220	zt(i) = zt(i) - zqbs(i) * (RLVTT+RLMLT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
221	ENDIF
222	zqbs(i)=0.
223	ENDIF
224
225
226	ENDDO ! loop on ngrid
227
228
229
230
231	! Outputs:
232	DO i = 1, ngrid
233	bsfl(i,k)=sedim(i)
234	dqbs_bs(i,k) = zqbs(i)-qbs(i,k)
235	dq_bs(i,k) = zq(i) - qv(i,k)
236	dtemp_bs(i,k) = zt(i) - temp(i,k)
237	ENDDO
238
239
240	ENDDO ! vertical loop
241
242
243	!surface bs flux
244	DO i = 1, ngrid
245	precip_bs(i) = sedim(i)
246	ENDDO
247
248
249	return
250
251	end subroutine blowing_snow_sublim_sedim
252	end module lmdz_blowing_snow_sublim_sedim

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format