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

Last change on this file since 4639 was 4523, checked in by evignon, 18 months ago

merge de la branche blowing snow vers la trunk
premiere tentative
Etienne
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.