source: LMDZ6/trunk/libf/phylmd/blowing_snow_sublim_sedim.F90 @ 4678

Last change on this file since 4678 was 4672, checked in by evignon, 9 months ago

inclusion des recents developpements sur la neige soufflee
File size: 7.7 KB
Line 
1subroutine blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,q,qbs,pplay,paprs,dtemp_bs,dq_bs,dqbs_bs,bsfl,precip_bs)
2
3!==============================================================================
4! Routine that calculates the evaporation and sedimentation of blowing snow
5! inspired by what is done in lscp_mod
6! Etienne Vignon, October 2022
7!==============================================================================
8
9
10use blowing_snow_ini_mod, only : coef_eva_bs,RTT,RD,RG,expo_eva_bs, fallv_bs
11use blowing_snow_ini_mod, only : RCPD, RLSTT, RLMLT, RLVTT, RVTMP2
12USE lmdz_lscp_tools, only : calc_qsat_ecmwf
13
14implicit none
15
16
17!++++++++++++++++++++++++++++++++++++++++++++++++++++
18! Declarations
19!++++++++++++++++++++++++++++++++++++++++++++++++++++
20
21!INPUT
22!=====
23
24integer, intent(in)                     :: ngrid,nlay
25real, intent(in)                        :: dtime
26real, intent(in), dimension(ngrid,nlay) :: temp
27real, intent(in), dimension(ngrid,nlay) :: q
28real, intent(in), dimension(ngrid,nlay) :: qbs
29real, intent(in), dimension(ngrid,nlay) :: pplay
30real, intent(in), dimension(ngrid,nlay+1) :: paprs
31
32
33
34! OUTPUT
35!========
36
37
38real, intent(out), dimension(ngrid,nlay) :: dtemp_bs
39real, intent(out), dimension(ngrid,nlay) :: dq_bs
40real, intent(out), dimension(ngrid,nlay) :: dqbs_bs
41real, intent(out), dimension(ngrid,nlay+1) :: bsfl
42real, intent(out), dimension(ngrid)      :: precip_bs
43
44
45! LOCAL
46!======
47
48
49integer                                  :: k,i,n
50real                                     :: zqev0, zqevi, zqevti, zcpair, zcpeau, dqbsmelt
51real                                     :: dqsedim,precbs
52real, dimension(ngrid)                   :: zt,zq,zqbs,qsi,dqsi,qsl, dqsl,qzero,sedim,sedimn
53real, dimension(ngrid)                   :: zqbsi,zmqc, zmair, zdz
54real, dimension(ngrid,nlay)              :: velo, zrho
55
56!++++++++++++++++++++++++++++++++++++++++++++++++++
57! Initialisation
58!++++++++++++++++++++++++++++++++++++++++++++++++++
59
60qzero(:)=0.
61dtemp_bs(:,:)=0.
62dq_bs(:,:)=0.
63dqbs_bs(:,:)=0.
64velo(:,:)=0.
65zt(:)=0.
66zq(:)=0.
67zqbs(:)=0.
68sedim(:)=0.
69
70! begin of top-down loop
71DO k = nlay, 1, -1
72   
73    DO i=1,ngrid
74        zt(i)=temp(i,k)
75        zq(i)=q(i,k)
76        zqbs(i)=qbs(i,k)
77    ENDDO
78
79
80    IF (k.LE.nlay-1) THEN
81
82        ! thermalization of blowing snow precip coming from above   
83        DO i = 1, ngrid
84            zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
85            ! RVTMP2=rcpv/rcpd-1
86            zcpair=RCPD*(1.0+RVTMP2*zq(i))
87            zcpeau=RCPD*RVTMP2
88            ! zmqc: precipitation mass that has to be thermalized with
89            ! layer's air so that precipitation at the ground has the
90            ! same temperature as the lowermost layer
91            zmqc(i) = (sedim(i))*dtime/zmair(i)
92            ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
93            zt(i) = ( (temp(i,k+1)+dtemp_bs(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) &
94             / (zcpair + zmqc(i)*zcpeau)
95        ENDDO
96    ELSE
97        DO i = 1, ngrid
98            zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
99            zmqc(i) = 0.
100        ENDDO
101
102    ENDIF
103
104
105
106    ! calulation saturation specific humidity
107    CALL CALC_QSAT_ECMWF(ngrid,zt(:),qzero(:),pplay(:,k),RTT,2,.false.,qsi(:),dqsi(:))
108    CALL CALC_QSAT_ECMWF(ngrid,zt(:),qzero(:),pplay(:,k),RTT,1,.false.,qsl(:),dqsl(:))
109    ! sublimation calculation
110    ! SUndqvist formula dP/dz=beta*(1-q/qsat)*sqrt(P)
111   
112    DO i = 1, ngrid
113       ! if sedimentation:
114       IF (sedim(i) .GT. 0.) THEN
115
116          IF (zt(i) .GT. RTT) THEN
117             ! if positive celcius temperature, we assume
118             ! that all the blowing snow melt and evaporate
119              zqevti=sedim(i)*(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG     
120             ! we ensure that the whole mesh does not exceed saturation wrt liquid
121              zqev0 = MAX(0.0, qsl(i)-zq(i))
122              zqevi = MIN(zqev0,zqevti)
123              ! New solid precipitation fluxes
124              sedimn(i) = Max(0.,sedim(i) - zqevi*(paprs(i,k)-paprs(i,k+1))/RG/dtime)
125              ! vapor, temperature, precip fluxes update
126              zq(i) = zq(i) - (sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime
127              zq(i) = max(0., zq(i))
128              !zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
129              !zqbs(i) = max(0., zqbs(i))
130              ! melting
131              zt(i) = zt(i)                             &
132                + (sedimn(i)-sedim(i))                      &
133                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
134                * RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
135              ! evaporation
136              zt(i) = zt(i)                             &
137                + (sedimn(i)-sedim(i))                      &
138                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
139                * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))     
140     
141     
142          ELSE     
143              zqevti = coef_eva_bs*(1.0-zq(i)/qsi(i))*(sedim(i)**expo_eva_bs) &
144                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
145              zqevti = MAX(0.0,MIN(zqevti,sedim(i)))*RG*dtime/(paprs(i,k)-paprs(i,k+1))
146
147              ! Sublimation limit: we ensure that the whole mesh does not exceed saturation wrt ice
148              zqev0 = MAX(0.0, qsi(i)-zq(i))
149              zqevi = MIN(zqev0,zqevti)
150
151              ! New solid precipitation fluxes
152              sedimn(i) = Max(0.,sedim(i) - zqevi*(paprs(i,k)-paprs(i,k+1))/RG/dtime)
153
154
155              ! vapor, temperature, precip fluxes update following sublimation
156              zq(i) = zq(i) - (sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime               
157              zq(i) = max(0., zq(i))
158              !zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
159              !zqbs(i) = max(0., zqbs(i))
160              zt(i) = zt(i)                             &
161                + (sedimn(i)-sedim(i))                      &
162                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
163                * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
164           ENDIF
165
166
167          ! sedim update
168          sedim(i)=sedimn(i)
169
170
171        ELSE
172           sedim(i)=0.
173        ENDIF ! if sedim > 0
174
175        zqbsi(i)=zqbs(i)
176
177    ENDDO ! loop on ngrid
178
179    ! Now sedimention scheme
180
181    ! exact resolution of the conservation equation for qbs with the updated flux from the top (after evap)
182    ! valid only if the fall velocity is constant
183
184    DO i = 1, ngrid
185   
186        zrho(i,k)  = pplay(i,k) / zt(i) / RD
187        zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
188        velo(i,k) = fallv_bs
189        zqbs(i) = zqbsi(i)*exp(-velo(i,k)/zdz(i)*dtime)+sedim(i)/zrho(i,k)/velo(i,k)
190        zqbs(i) = max(zqbs(i),0.)
191        ! flux update
192        sedim(i) = sedim(i) + zrho(i,k)*zdz(i)/dtime*(zqbsi(i)-zqbs(i))
193        sedim(i) = max(0.,sedim(i))
194 
195    ENDDO
196
197! old version with bugs
198!    DO n = 1, niter_bs
199!       DO i = 1, ngrid
200!         zrho(i,k)  = pplay(i,k) / zt(i) / RD
201!         zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
202!         velo(i,k) = fallv_bs
203!         dqsedim = dtime/REAL(niter_bs)/zdz(i)*zqbs(i)*velo(i,k)   ! dqice/dt=1/rho*d(rho*wice*qice)/dz
204!         precbs   = MIN(MAX(dqsedim,0.0),zqbs(i))
205!         zqbs(i) = MAX(zqbs(i)-1.*precbs  , 0.0)
206!       ENDDO !loop on ngrid
207!    ENDDO ! loop on niter_bs
208!    ! add to non sublimated precip
209!    DO i=1,ngrid
210!       sedim(i) = sedim(i)+max(0.,zqbsi(i)-zqbs(i))*(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
211!    ENDDO
212!
213
214
215
216    ! Outputs:
217    DO i = 1, ngrid
218        bsfl(i,k)=sedim(i)
219        dqbs_bs(i,k) = zqbs(i)-qbs(i,k)
220        dq_bs(i,k) = zq(i) - q(i,k)
221        dtemp_bs(i,k) = zt(i) - temp(i,k)
222    ENDDO
223
224
225ENDDO ! vertical loop
226
227
228!surface bs flux
229DO i = 1, ngrid
230  precip_bs(i) = sedim(i)
231ENDDO
232
233
234return
235
236end subroutine blowing_snow_sublim_sedim
Note: See TracBrowser for help on using the repository browser.