source: LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_sublim_sedim.f90 @ 5912

Last change on this file since 5912 was 5912, checked in by evignon, 8 hours ago

changes in blowing snow routines to comply with the revised version of the paper
File size: 10.7 KB
Line 
1module lmdz_blowing_snow_sublim_sedim
2
3contains
4
5!==============================================================================       
6subroutine blowing_snow_sublim_sedim(ngrid,nlay,dtime,ustar,temp,qv,qb,pplay,paprs,dtemp_bs,dqv_bs,dqb_bs,bsfl,precip_bs)
7
8!==============================================================================
9! Routine that calculates the sublimation, melting and sedimentation
10! of blowing snow
11! Reference: Vignon et al 2025, GMD, https://doi.org/10.5194/egusphere-2025-2871
12!
13! Etienne Vignon, October 2022
14!==============================================================================
15
16
17use lmdz_blowing_snow_ini, only : iflag_sublim_bs, iflag_sedim_bs, coef_sub_bs,RTT,RD,RG,fallv_bs
18use lmdz_blowing_snow_ini, only : qbmin, RCPD, RLSTT, RLMLT, RLVTT, RVTMP2, RV, RPI
19use lmdz_blowing_snow_ini, only : tbsmelt, taumeltbs0, rhobs, r_bs
20USE lmdz_lscp_tools, only : calc_qsat_ecmwf
21
22implicit none
23
24
25!++++++++++++++++++++++++++++++++++++++++++++++++++++
26! Declarations
27!++++++++++++++++++++++++++++++++++++++++++++++++++++
28
29!INPUT
30!=====
31
32integer, intent(in)                     :: ngrid ! number of horizontal grid points
33integer, intent(in)                     :: nlay  ! number of vertical layers
34real, intent(in)                        :: dtime ! physics time step [s]
35real, intent(in), dimension(ngrid)      :: ustar ! surface friction velocity [m/s]
36real, intent(in), dimension(ngrid,nlay) :: temp  ! temperature of the air [K]
37real, intent(in), dimension(ngrid,nlay) :: qv    ! specific content of water [kg/kg]
38real, intent(in), dimension(ngrid,nlay) :: qb    ! specific content of blowing snow [kg/kg]
39real, intent(in), dimension(ngrid,nlay) :: pplay ! pressure at the middle of layers [Pa]
40real, intent(in), dimension(ngrid,nlay+1) :: paprs ! pressure at the layer bottom interface [Pa]
41
42
43
44! OUTPUT
45!========
46
47
48real, intent(out), dimension(ngrid,nlay) :: dtemp_bs ! temperature tendency [K/s]
49real, intent(out), dimension(ngrid,nlay) :: dqv_bs   ! water vapor tendency [kg/kg/s]
50real, intent(out), dimension(ngrid,nlay) :: dqb_bs   ! blowing snow mass tendancy [kg/kg/s]
51real, intent(out), dimension(ngrid,nlay+1) :: bsfl   ! vertical profile of blowing snow vertical flux [kg/m2/s]
52real, intent(out), dimension(ngrid)      :: precip_bs ! surface sedimentation flux of blowing snow [kg/s/m2]
53
54! LOCAL
55!======
56
57
58integer                                  :: k,i,n
59real                                     :: cpd, cpw, dqbsub, maxdqbsub, dqbmelt, zmair
60real                                     :: dqsedim,precbs, dqvmelt, zmelt, taumeltbs
61real                                     :: maxdqvmelt, rhoair, dz, dqbsedim
62real                                     :: delta_p, b_p, a_p, c_p, c_sub, qvsub
63real                                     :: p0, T0, Dv, Aprime, Bprime, Ka, radius0, radius
64real, dimension(ngrid)                   :: ztemp,zqv,zqb,zpres,qsi,dqsi,qsl,dqsl,qzero,sedim
65real, dimension(ngrid)                   :: zpaprsup, zpaprsdn, ztemp_up, zqb_up, zvelo
66real, dimension(ngrid,nlay)              :: temp_seri, qb_seri, qv_seri, zz
67
68!++++++++++++++++++++++++++++++++++++++++++++++++++
69! Initialisation
70!++++++++++++++++++++++++++++++++++++++++++++++++++
71
72qzero(:)=0.
73dtemp_bs(:,:)=0.
74dqv_bs(:,:)=0.
75dqb_bs(:,:)=0.
76zvelo(:)=0.
77sedim(:)=0.
78precip_bs(:)=0.
79bsfl(:,:)=0.
80zz(:,:)=0.
81
82
83p0=101325.0    ! ref pressure
84
85
86DO k=1,nlay
87   DO i=1,ngrid
88      temp_seri(i,k)=temp(i,k)
89      qv_seri(i,k)=qv(i,k)
90      qb_seri(i,k)=qb(i,k)
91      zz(i,k)=zz(i,k)+(paprs(i,k)-paprs(i,k+1)) / (pplay(i,k)/RD/temp(i,k)*RG)
92   ENDDO
93ENDDO
94
95
96! Sedimentation scheme
97!----------------------
98
99IF (iflag_sedim_bs .GT. 0) THEN
100! begin of top-down loop
101DO k = nlay, 1, -1
102   
103    DO i=1,ngrid
104        ztemp(i)=temp_seri(i,k)
105        zqv(i)=qv_seri(i,k)
106        zqb(i)=qb_seri(i,k)
107        zpres(i)=pplay(i,k)
108        zpaprsup(i)=paprs(i,k+1)
109        zpaprsdn(i)=paprs(i,k)
110    ENDDO
111
112    ! thermalization of blowing snow precip coming from above     
113    IF (k.LE.nlay-1) THEN
114
115        DO i = 1, ngrid
116            zmair=(zpaprsdn(i)-zpaprsup(i))/RG
117            ! RVTMP2=rcpv/rcpd-1
118            cpd=RCPD*(1.0+RVTMP2*zqv(i))
119            cpw=RCPD*RVTMP2
120            ! zqb_up: blowing snow mass that has to be thermalized with
121            ! layer's air so that precipitation at the ground has the
122            ! same temperature as the lowermost layer
123            zqb_up(i) = sedim(i)*dtime/zmair
124            ztemp_up(i)=temp_seri(i,k+1)
125
126            ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
127            ztemp(i) = ( ztemp_up(i)*zqb_up(i)*cpw + cpd*ztemp(i) ) &
128                  / (cpd + zqb_up(i)*cpw)
129        ENDDO
130
131     ENDIF
132
133     DO i = 1, ngrid
134
135       rhoair  = zpres(i) / ztemp(i) / RD
136       dz      = (zpaprsdn(i)-zpaprsup(i)) / (rhoair*RG)
137       ! BS fall velocity assumed to be constant for now
138       zvelo(i) = fallv_bs
139       ! dqb/dt_sedim=1/rho(sedim/dz - rho w qb/dz)
140       ! implicit resolution
141       dqbsedim = (sedim(i)/rhoair/dz*dtime+zqb(i))/(1.+zvelo(i)*dtime/dz) - zqb(i)
142       ! flux and dqb update
143       zqb(i)=zqb(i)+dqbsedim
144       !sedim(i) = sedim(i)-dqbsedim/dtime*rhoair*dz
145       sedim(i)=rhoair*zvelo(i)*zqb(i)
146
147
148       ! variables update:
149       bsfl(i,k)=sedim(i)
150
151       qb_seri(i,k) = zqb(i)
152       qv_seri(i,k) = zqv(i)
153       temp_seri(i,k) = ztemp(i)
154
155     ENDDO ! Loop on ngrid
156
157
158ENDDO ! vertical loop
159
160
161!surface bs flux
162DO i = 1, ngrid
163  precip_bs(i) = sedim(i)
164ENDDO
165
166ENDIF
167
168
169
170
171!+++++++++++++++++++++++++++++++++++++++++++++++
172! Sublimation and melting
173!++++++++++++++++++++++++++++++++++++++++++++++
174IF (iflag_sublim_bs .GT. 0) THEN
175
176DO k = 1, nlay
177
178    DO i=1,ngrid
179        ztemp(i)=temp_seri(i,k)
180        zqv(i)=qv_seri(i,k)
181        zqb(i)=qb_seri(i,k)
182        zpres(i)=pplay(i,k)
183        zpaprsup(i)=paprs(i,k+1)
184        zpaprsdn(i)=paprs(i,k)
185
186    ! computation of blowing snow mean radius. Either constant value = r_bs if >0
187    ! or use of the height-dependent formulation from Sauter et al. 2013 and Saigger et al. 2024
188        IF (r_bs .GT. 0.) THEN
189            radius=r_bs
190        ELSE
191            radius0=0.5*(ustar(i)*(7.8e-6)/0.036+31.e-6)
192            radius=radius0*zz(i,k)**(-0.258)
193        ENDIF
194
195    ENDDO
196
197   ! calulation saturation specific humidity
198    CALL CALC_QSAT_ECMWF(ngrid,ztemp(:),qzero(:),zpres(:),RTT,2,.false.,qsi(:),dqsi(:))
199
200
201    DO i = 1, ngrid
202
203          rhoair  = zpres(i) / ztemp(i) / RD
204          dz      = (zpaprsdn(i)-zpaprsup(i)) / (rhoair*RG)
205          ! BS fall velocity assumed to be constant for now
206          zvelo(i) = fallv_bs
207
208
209          IF (ztemp(i) .GT. RTT) THEN
210
211             ! if temperature is positive, we assume that part of the blowing snow
212             ! already present  melts and evaporates with a typical time
213             ! constant taumeltbs
214
215             taumeltbs=taumeltbs0*exp(-max(0.,(ztemp(i)-RTT)/(tbsmelt-RTT)))
216             dqvmelt=min(zqb(i),-1.*zqb(i)*(exp(-dtime/taumeltbs)-1.))
217             maxdqvmelt= max(RCPD*(1.0+RVTMP2*(zqv(i)+zqb(i)))*(ztemp(i)-RTT)/(RLMLT+RLVTT),0.)
218             dqvmelt=min(dqvmelt,maxdqvmelt)
219             ! qv update, melting + evaporation
220             zqv(i) = zqv(i) + dqvmelt
221             ! temp update melting + evaporation
222             ztemp(i) = ztemp(i) - dqvmelt * (RLMLT+RLVTT)/RCPD/(1.0+RVTMP2*(zqv(i)+zqb(i)))
223             ! qb update melting + evaporation
224             zqb(i)=zqb(i)-dqvmelt
225
226          ELSE
227              ! negative celcius temperature     
228              ! Sublimation scheme   
229
230
231              ! Sublimation formulation for ice crystals from Pruppacher & Klett, Rutledge & Hobbs 1983
232              ! assuming monodispered crystal distrib
233              ! dqb/dt_sub=-coef_sub_bs*(1-qv/qsi)*nc*8*radius/(Aprime+Bprime)
234              ! assuming Mi=rhobs*4/3*pi*radius**3
235              ! qb=nc*Mi -> nc=qb/Mi
236              ! dqb/dt_sub=-coef_sub_bs*(1-qv/qsi)*6*qb/(rhobs*pi*radius**2)/(Aprime+Bprime)
237              ! dqb/dt_sub=-c_sub(1-qv/qsi)*qb
238              ! c_sub=coef_sub_bs*6/(rhobs*pi*radius**2)/(Aprime+Bprime)
239              !
240              ! Note the strong coupling between specific contents of water vapor and blowing snow during sublimation
241              ! equations dqv/dt_sub and dqb/dt_sub must be solved jointly to prevent irrealistic increase of water vapor
242              ! at typical physics time steps
243              ! we thus solve the differential equation using an implicit scheme for both qb and qv
244 
245              ! we do not consider deposition, only sublimation
246              IF (zqv(i) .LT. qsi(i)) THEN
247                 rhoair=zpres(i)/ztemp(i)/RD
248                 Dv=0.0001*0.211*(p0/zpres(i))*((ztemp(i)/RTT)**1.94)           ! water vapor diffusivity in air, SI
249                 Ka=(5.69+0.017*(ztemp(i)-RTT))*1.e-5*100.*4.184                ! thermal conductivity of the air, SI
250                 Aprime=RLSTT/Ka/ztemp(i)*(RLSTT/RV/ztemp(i) -1.)
251                 Bprime=1./(rhoair*Dv*qsi(i))
252                 c_sub=coef_sub_bs*3./(rhobs*radius**2)/(Aprime+Bprime)
253                 c_p=-zqb(i)
254                 b_p=1.+c_sub*dtime-c_sub*dtime/qsi(i)*zqb(i)-c_sub*dtime/qsi(i)*zqv(i)
255                 a_p=c_sub*dtime/qsi(i)
256                 delta_p=(b_p**2)-4.*a_p*c_p 
257                 dqbsub=(-b_p+sqrt(delta_p))/(2.*a_p) - zqb(i)       
258                 dqbsub = MIN(0.0,MAX(dqbsub,-zqb(i)))
259                 ! Sublimation limit: we ensure that the whole mesh does not exceed saturation wrt ice
260                 maxdqbsub = MAX(0.0, qsi(i)-zqv(i))
261                 dqbsub = MAX(dqbsub,-maxdqbsub)
262              ELSE
263                 dqbsub=0.
264              ENDIF
265
266              ! vapor, temperature, precip fluxes update following sublimation
267              zqv(i) = zqv(i) - dqbsub
268              zqb(i) = zqb(i) + dqbsub
269              ztemp(i) = ztemp(i) + dqbsub*RLSTT/RCPD/(1.0+RVTMP2*(zqv(i)+zqb(i)))
270
271          ENDIF
272
273          ! if qb<qbmin, sublimate or melt and evaporate qb
274          ! see Gerber et al. 2023, JGR Atmos for the choice of qbmin
275
276          IF (zqb(i) .LT. qbmin) THEN
277              zqv(i) = zqv(i)+zqb(i)
278              IF (ztemp(i) .LT. RTT) THEN
279                 ztemp(i) = ztemp(i) - zqb(i) * RLSTT/RCPD/(1.0+RVTMP2*(zqv(i)))
280              ELSE
281                 ztemp(i) = ztemp(i) - zqb(i) * (RLVTT+RLMLT)/RCPD/(1.0+RVTMP2*(zqv(i)))
282              ENDIF
283              zqb(i)=0.
284          ENDIF
285
286     ! variables update
287     temp_seri(i,k)=ztemp(i)
288     qv_seri(i,k)=zqv(i)
289     qb_seri(i,k)=zqb(i)
290     ENDDO
291ENDDO
292
293ENDIF
294
295
296! OUTPUTS
297!++++++++++
298
299! 3D variables
300DO k=1,nlay
301   DO i=1,ngrid
302        dqb_bs(i,k) = qb_seri(i,k) - qb(i,k)
303        dqv_bs(i,k) = qv_seri(i,k) - qv(i,k)
304        dtemp_bs(i,k) = temp_seri(i,k) - temp(i,k)
305   ENDDO
306ENDDO
307
308
309
310return
311
312end subroutine blowing_snow_sublim_sedim
313!==============================================================================
314
315
316
317end module lmdz_blowing_snow_sublim_sedim
Note: See TracBrowser for help on using the repository browser.