source: LMDZ6/trunk/libf/phylmd/lmdz_blowing_snow_sublim_sedim.F90 @ 4821

Last change on this file since 4821 was 4821, checked in by evignon, 3 months ago

reecriture de la routine de sublimation de neige soufflee en vue
du debut d'un stage de M2 sur le sujet prochainement
File size: 6.3 KB
Line 
1module lmdz_blowing_snow_sublim_sedim
2
3contains
4subroutine blowing_snow_sublim_sedim(ngrid,nlay,dtime,temp,qv,qbs,pplay,paprs,dtemp_bs,dq_bs,dqbs_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 : coef_eva_bs,RTT,RD,RG,expo_eva_bs, fallv_bs, qbsmin
14use lmdz_blowing_snow_ini, only : RCPD, RLSTT, RLMLT, RLVTT, RVTMP2, tbsmelt, taumeltbs0
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) :: qbs
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) :: dq_bs
43real, intent(out), dimension(ngrid,nlay) :: dqbs_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, dqsub, maxdqsub, dqbsmelt
54real                                     :: dqsedim,precbs, dqmelt, zmelt, taumeltbs
55real                                     :: maxdqmelt, rhoair, dz
56real, dimension(ngrid)                   :: zt,zq,zqbs,qsi,dqsi,qsl, dqsl,qzero,sedim,sedimn
57real, dimension(ngrid)                   :: zqbsi, zqbs_up, zmair
58real, dimension(ngrid)                   :: zvelo
59
60!++++++++++++++++++++++++++++++++++++++++++++++++++
61! Initialisation
62!++++++++++++++++++++++++++++++++++++++++++++++++++
63
64qzero(:)=0.
65dtemp_bs(:,:)=0.
66dq_bs(:,:)=0.
67dqbs_bs(:,:)=0.
68zvelo(:)=0.
69zt(:)=0.
70zq(:)=0.
71zqbs(:)=0.
72sedim(:)=0.
73sedimn(:)=0.
74
75! begin of top-down loop
76DO k = nlay, 1, -1
77   
78    DO i=1,ngrid
79        zt(i)=temp(i,k)
80        zq(i)=qv(i,k)
81        zqbs(i)=qbs(i,k)
82    ENDDO
83
84     ! thermalization of blowing snow precip coming from above     
85    IF (k.LE.nlay-1) THEN
86
87        DO i = 1, ngrid
88            zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG
89            ! RVTMP2=rcpv/rcpd-1
90            cpd=RCPD*(1.0+RVTMP2*zq(i))
91            cpw=RCPD*RVTMP2
92            ! zqbs_up: blowing snow mass that has to be thermalized with
93            ! layer's air so that precipitation at the ground has the
94            ! same temperature as the lowermost layer
95            zqbs_up(i) = sedim(i)*dtime/zmair(i)
96            ! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer
97            zt(i) = ( (temp(i,k+1)+dtemp_bs(i,k+1))*zqbs_up(i)*cpw + cpd*zt(i) ) &
98             / (cpd + zqbs_up(i)*cpw)
99        ENDDO
100    ENDIF
101
102
103
104    ! calulation saturation specific humidity
105    CALL CALC_QSAT_ECMWF(ngrid,zt(:),qzero(:),pplay(:,k),RTT,2,.false.,qsi(:),dqsi(:))
106    CALL CALC_QSAT_ECMWF(ngrid,zt(:),qzero(:),pplay(:,k),RTT,1,.false.,qsl(:),dqsl(:))
107   
108   
109    ! 3 processes: melting, sublimation and precipitation of blowing snow
110    DO i = 1, ngrid
111
112          IF (zt(i) .GT. RTT) THEN
113             
114             ! if temperature is positive, we assume that part of the blowing snow
115             ! already present  melts and evaporates with a typical time
116             ! constant taumeltbs
117           
118             taumeltbs=taumeltbs0*exp(-max(0.,(zt(i)-RTT)/(tbsmelt-RTT)))
119             dqmelt=min(zqbs(i),zqbs(i)/taumeltbs*dtime)
120             maxdqmelt= max(RCPD*(1.0+RVTMP2*(zq(i)+zqbs(i)))*(zt(i)-RTT)/(RLMLT+RLVTT),0.)
121             dqmelt=min(max(dqmelt,0.),maxdqmelt)
122             ! q update, melting + evaporation
123             zq(i) = zq(i) + dqmelt
124             ! temp update melting + evaporation
125             zt(i) = zt(i) - dqmelt * (RLMLT+RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
126             ! qbs update melting + evaporation
127             zqbs(i)=zqbs(i)-dqmelt
128
129          ELSE
130              ! negative celcius temperature     
131              ! Sublimation scheme   
132
133              rhoair=pplay(i,k)/zt(i)/RD
134              dqsub = 1./rhoair*coef_eva_bs*(1.0-zq(i)/qsi(i))*((rhoair*zvelo(i)*zqbs(i))**expo_eva_bs)*dtime
135              dqsub = MAX(0.0,MIN(dqsub,zqbs(i)))
136
137              ! Sublimation limit: we ensure that the whole mesh does not exceed saturation wrt ice
138              maxdqsub = MAX(0.0, qsi(i)-zq(i))
139              dqsub = MIN(dqsub,maxdqsub)
140
141
142              ! vapor, temperature, precip fluxes update following sublimation
143              zq(i) = zq(i) + dqsub     
144              zqbs(i)=zqbs(i)-dqsub           
145              zt(i) = zt(i) - dqsub*RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
146
147          ENDIF
148
149        ! Sedimentation scheme
150           
151        rhoair  = pplay(i,k) / zt(i) / RD
152        dz      = (paprs(i,k)-paprs(i,k+1)) / (rhoair*RG)
153        ! BS fall velocity assumed to be constant for now
154        zvelo(i) = fallv_bs
155
156        sedimn(i) = rhoair*zqbs(i)*zvelo(i)
157
158        ! exact numerical resolution of sedimentation
159        ! assuming fall velocity is constant
160
161        zqbs(i) = zqbs(i)*exp(-zvelo(i)/dz*dtime)+sedim(i)/rhoair/zvelo(i)
162
163        ! flux update
164        sedim(i) = sedimn(i)
165
166
167
168        ! if qbs<qbsmin, sublimate or melt and evaporate qbs
169        ! see Gerber et al. 2023, JGR Atmos for the choice of qbsmin
170        IF (zqbs(i) .LT. qbsmin) THEN                       
171              zq(i) = zq(i)+zqbs(i)
172              IF (zt(i) .LT. RTT) THEN
173                 zt(i) = zt(i) - zqbs(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
174              ELSE
175               zt(i) = zt(i) - zqbs(i) * (RLVTT+RLMLT)/RCPD/(1.0+RVTMP2*(zq(i)+zqbs(i)))
176              ENDIF
177              zqbs(i)=0.
178        ENDIF
179
180
181    ! Outputs:
182        bsfl(i,k)=sedim(i)
183        dqbs_bs(i,k) = zqbs(i)-qbs(i,k)
184        dq_bs(i,k) = zq(i) - qv(i,k)
185        dtemp_bs(i,k) = zt(i) - temp(i,k)
186
187    ENDDO ! Loop on ngrid
188
189
190ENDDO ! vertical loop
191
192
193!surface bs flux
194DO i = 1, ngrid
195  precip_bs(i) = sedim(i)
196ENDDO
197
198
199return
200
201end subroutine blowing_snow_sublim_sedim
202end module lmdz_blowing_snow_sublim_sedim
Note: See TracBrowser for help on using the repository browser.