source: LMDZ6/branches/contrails/libf/phylmd/lmdz_blowing_snow_sublim_sedim.f90 @ 5450

Last change on this file since 5450 was 5268, checked in by abarral, 2 months ago

.f90 <-> .F90 depending on cpp key use
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.