source: LMDZ6/branches/lmdz-snow/libf/phylmd/lmdz_blowing_snow_sublim_sedim.f90 @ 5872

Last change on this file since 5872 was 5872, checked in by evignon, 8 days ago

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