source: LMDZ6/branches/blowing_snow/libf/phylmd/blowing_snow_sublim_sedim.F90

Last change on this file was 4485, checked in by evignon, 21 months ago

premier commit pour l'ajout de la neige soufflee sur la nouvelle branche
File size: 6.9 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, niter_bs
12USE lscp_tools_mod, 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              zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
128              ! melting
129              zt(i) = zt(i)                             &
130                + (sedimn(i)-sedim(i))                      &
131                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
132                * RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
133              ! evaporation
134              zt(i) = zt(i)                             &
135                + (sedimn(i)-sedim(i))                      &
136                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
137                * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))     
138     
139     
140          ELSE     
141              zqevti = coef_eva_bs*(1.0-zq(i)/qsi(i))*(sedim(i)**expo_eva_bs) &
142                    *(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG
143              zqevti = MAX(0.0,MIN(zqevti,sedim(i)))*RG*dtime/(paprs(i,k)-paprs(i,k+1))
144
145              ! Sublimation limit: we ensure that the whole mesh does not exceed saturation wrt ice
146              zqev0 = MAX(0.0, qsi(i)-zq(i))
147              zqevi = MIN(zqev0,zqevti)
148
149              ! New solid precipitation fluxes
150              sedimn(i) = Max(0.,sedim(i) - zqevi*(paprs(i,k)-paprs(i,k+1))/RG/dtime)
151
152
153              ! vapor, temperature, precip fluxes update
154              zq(i) = zq(i) - (sedimn(i)-sedim(i))*(RG/(paprs(i,k)-paprs(i,k+1)))*dtime               
155              zqbs(i) = zqbs(i) + (sedimn(i)-sedim(i)) * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime
156              zt(i) = zt(i)                             &
157                + (sedimn(i)-sedim(i))                      &
158                * (RG/(paprs(i,k)-paprs(i,k+1)))*dtime    &
159                * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
160           ENDIF
161
162
163          ! sedim update
164          sedim(i)=sedimn(i)
165
166
167        ELSE
168           sedim(i)=0.
169        ENDIF ! if sedim > 0
170
171        zqbsi(i)=zqbs(i)
172
173    ENDDO ! loop on ngrid
174
175    ! Now sedimention scheme (alike that in lscp)
176    DO n = 1, niter_bs
177       DO i = 1, ngrid
178         zrho(i,k)  = pplay(i,k) / zt(i) / RD
179         zdz(i)   = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG)
180         velo(i,k) = fallv_bs
181         dqsedim = dtime/REAL(niter_bs)/zdz(i)*zqbs(i)*velo(i,k)   ! dqice/dt=1/rho*d(rho*wice*qice)/dz
182         precbs   = MIN(MAX(dqsedim,0.0),zqbs(i))
183         zqbs(i) = MAX(zqbs(i)-1.*precbs  , 0.0)
184       ENDDO !loop on ngrid
185    ENDDO ! loop on niter_bs
186
187    ! add to non sublimated precip
188    DO i=1,ngrid
189    sedim(i) = sedim(i)+max(0.,zqbsi(i)-zqbs(i))*(paprs(i,k)-paprs(i,k+1))/(RG*dtime)
190    ENDDO
191
192    ! Outputs:
193    DO i = 1, ngrid
194        bsfl(i,k)=sedim(i)
195        dqbs_bs(i,k) = zqbs(i)-qbs(i,k)
196        dq_bs(i,k) = zq(i) - q(i,k)
197        dtemp_bs(i,k) = zt(i) - temp(i,k)
198    ENDDO
199
200
201ENDDO ! vertical loop
202
203
204!surface bs flux
205DO i = 1, ngrid
206  precip_bs(i) = sedim(i)
207ENDDO
208
209
210return
211
212end subroutine blowing_snow_sublim_sedim
Note: See TracBrowser for help on using the repository browser.