source: LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_sublim_sedim.F90 @ 4835

Last change on this file since 4835 was 4835, checked in by evignon, 3 months ago

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