source: LMDZ6/branches/cirrus/libf/phylmd/lmdz_blowing_snow_sublim_sedim.F90 @ 5507

Last change on this file since 5507 was 4916, checked in by evignon, 10 months ago

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