1 |
|
|
MODULE calcratqs_multi_mod |
2 |
|
|
|
3 |
|
|
!============================================= |
4 |
|
|
! module containing subroutines that take |
5 |
|
|
! into account the effect of convection, orography, |
6 |
|
|
! surface heterogeneities and subgrid-scale |
7 |
|
|
! turbulence on ratqs, i.e. on the width of the |
8 |
|
|
! total water subgrid distribution. |
9 |
|
|
!============================================= |
10 |
|
|
|
11 |
|
|
IMPLICIT NONE |
12 |
|
|
|
13 |
|
|
! Include |
14 |
|
|
!============================================= |
15 |
|
|
INCLUDE "YOETHF.h" |
16 |
|
|
INCLUDE "YOMCST.h" |
17 |
|
|
|
18 |
|
|
|
19 |
|
|
CONTAINS |
20 |
|
|
|
21 |
|
|
|
22 |
|
|
!======================================================================== |
23 |
|
|
SUBROUTINE calcratqs_inter(klon,klev,iflag_ratqs,pdtphys, & |
24 |
|
|
ratqsbas, wake_deltaq, wake_s, q_seri,qtc_cv, sigt_cv, & |
25 |
|
|
ratqs_inter) |
26 |
|
|
USE ioipsl_getin_p_mod, ONLY : getin_p |
27 |
|
|
implicit none |
28 |
|
|
|
29 |
|
|
!======================================================================== |
30 |
|
|
! L. d'Alençon, 25/02/2021 |
31 |
|
|
! Cette subroutine calcule une valeur de ratqsbas interactive dépendant de la présence de poches froides dans l'environnement. |
32 |
|
|
! Elle est appelée par la subroutine calcratqs lorsque iflag_ratqs = 10. |
33 |
|
|
!======================================================================== |
34 |
|
|
|
35 |
|
|
! Declarations |
36 |
|
|
|
37 |
|
|
|
38 |
|
|
LOGICAL, SAVE :: first = .TRUE. |
39 |
|
|
!$OMP THREADPRIVATE(first) |
40 |
|
|
! Input |
41 |
|
|
integer,intent(in) :: klon,klev,iflag_ratqs |
42 |
|
|
real,intent(in) :: pdtphys,ratqsbas |
43 |
|
|
real, dimension(klon,klev),intent(in) :: wake_deltaq, q_seri,qtc_cv, sigt_cv |
44 |
|
|
real, dimension(klon),intent(in) :: wake_s |
45 |
|
|
|
46 |
|
|
! Output |
47 |
|
|
real, dimension(klon,klev),intent(inout) :: ratqs_inter |
48 |
|
|
|
49 |
|
|
! local |
50 |
|
|
integer i,k |
51 |
|
|
real, dimension(klon,klev) :: wake_dq |
52 |
|
|
REAL, SAVE :: a_ratqs_cv |
53 |
|
|
!$OMP THREADPRIVATE(a_ratqs_cv) |
54 |
|
|
REAL, SAVE :: tau_ratqs_wake |
55 |
|
|
!$OMP THREADPRIVATE(tau_ratqs_wake) |
56 |
|
|
REAL, SAVE :: a_ratqs_wake |
57 |
|
|
!$OMP THREADPRIVATE(a_ratqs_wake) |
58 |
|
|
real, dimension(klon) :: max_wake_dq, max_dqconv,max_sigt |
59 |
|
|
!------------------------------------------------------------------------- |
60 |
|
|
! Caclul de ratqs_inter |
61 |
|
|
!------------------------------------------------------------------------- |
62 |
|
|
|
63 |
|
|
! |
64 |
|
|
if (first) then |
65 |
|
|
tau_ratqs_wake = 3600. ! temps de relaxation de la variabilité |
66 |
|
|
a_ratqs_wake = 3. ! paramètre pilotant l'importance du terme dépendant des poches froides |
67 |
|
|
a_ratqs_cv = 0.5 |
68 |
|
|
CALL getin_p('tau_ratqs_wake', tau_ratqs_wake) |
69 |
|
|
CALL getin_p('a_ratqs_wake', a_ratqs_wake) |
70 |
|
|
CALL getin_p('a_ratqs_cv', a_ratqs_cv) |
71 |
|
|
first=.false. |
72 |
|
|
endif |
73 |
|
|
max_wake_dq(:) = 0. |
74 |
|
|
max_dqconv (:) = 0 |
75 |
|
|
max_sigt(:) = 0. |
76 |
|
|
if (iflag_ratqs.eq.10) then |
77 |
|
|
do k=1,klev |
78 |
|
|
do i=1,klon |
79 |
|
|
max_wake_dq(i) = max(abs(wake_deltaq(i,k)),max_wake_dq(i)) |
80 |
|
|
max_sigt(i) = max(abs(sigt_cv(i,k)),max_sigt(i)) |
81 |
|
|
max_dqconv(i) = max(abs(q_seri(i,k) - qtc_cv(i,k)),max_dqconv(i)) |
82 |
|
|
enddo |
83 |
|
|
enddo |
84 |
|
|
|
85 |
|
|
do k=1,klev |
86 |
|
|
do i=1,klon |
87 |
|
|
ratqs_inter(i,k)= ratqs_inter(i,k)*exp(-pdtphys/tau_ratqs_wake) + & |
88 |
|
|
a_ratqs_wake*(max_wake_dq(i)*(wake_s(i)**0.5/(1.-wake_s(i))))*(1.-exp(-pdtphys/tau_ratqs_wake))/q_seri(i,1) |
89 |
|
|
if (ratqs_inter(i,k)<ratqsbas) then |
90 |
|
|
ratqs_inter(i,k) = ratqsbas |
91 |
|
|
endif |
92 |
|
|
enddo |
93 |
|
|
enddo |
94 |
|
|
endif |
95 |
|
|
|
96 |
|
|
if (iflag_ratqs.eq.11) then |
97 |
|
|
do k=1,klev |
98 |
|
|
do i=1,klon |
99 |
|
|
max_wake_dq(i) = max(abs(wake_deltaq(i,k)),max_wake_dq(i)) |
100 |
|
|
max_sigt(i) = max(abs(sigt_cv(i,k)),max_sigt(i)) |
101 |
|
|
max_dqconv(i) = max(abs(q_seri(i,k) - qtc_cv(i,k)),max_dqconv(i)) |
102 |
|
|
enddo |
103 |
|
|
enddo |
104 |
|
|
|
105 |
|
|
do k=1,klev |
106 |
|
|
do i=1,klon |
107 |
|
|
ratqs_inter(i,k)= ratqs_inter(i,k)*exp(-pdtphys/tau_ratqs_wake) + & |
108 |
|
|
a_ratqs_wake*(max_wake_dq(i)*(wake_s(i)**0.5/(1.-wake_s(i))))*(1.-exp(-pdtphys/tau_ratqs_wake))/q_seri(i,1) + & |
109 |
|
|
a_ratqs_cv*max_dqconv(i)*max_sigt(i)*(1.-exp(-pdtphys/tau_ratqs_wake))/q_seri(i,1) |
110 |
|
|
! if (ratqs_inter(i,k)>0) then |
111 |
|
|
! ratqs_inter(i,k) = abs(q_seri(i,k) - qtc_cv(i,k)) |
112 |
|
|
! endif |
113 |
|
|
enddo |
114 |
|
|
enddo |
115 |
|
|
endif |
116 |
|
|
return |
117 |
|
|
end |
118 |
|
|
|
119 |
|
|
|
120 |
|
|
!------------------------------------------------------------------ |
121 |
|
|
SUBROUTINE calcratqs_oro(klon,klev,qsat,temp,pplay,paprs,ratqs_oro) |
122 |
|
|
|
123 |
|
|
! Etienne Vignon, November 2021: effect of subgrid orography on ratqs |
124 |
|
|
|
125 |
|
|
USE phys_state_var_mod, ONLY: zstd |
126 |
|
|
USE phys_state_var_mod, ONLY: pctsrf |
127 |
|
|
USE indice_sol_mod, only: nbsrf, is_lic, is_ter |
128 |
|
|
|
129 |
|
|
IMPLICIT NONE |
130 |
|
|
|
131 |
|
|
! Declarations |
132 |
|
|
!-------------- |
133 |
|
|
|
134 |
|
|
! INPUTS |
135 |
|
|
|
136 |
|
|
INTEGER, INTENT(IN) :: klon ! number of horizontal grid points |
137 |
|
|
INTEGER, INTENT(IN) :: klev ! number of vertical layers |
138 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: qsat ! saturation specific humidity [kg/kg] |
139 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! air temperature [K] |
140 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! air pressure, layer's center [Pa] |
141 |
|
|
REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! air pressure, lower inteface [Pa] |
142 |
|
|
|
143 |
|
|
! OUTPUTS |
144 |
|
|
|
145 |
|
|
REAL, DIMENSION(klon,klev), INTENT(out) :: ratqs_oro ! ratqs profile due to subgrid orography |
146 |
|
|
|
147 |
|
|
|
148 |
|
|
! LOCAL |
149 |
|
|
|
150 |
|
|
INTEGER :: i,k |
151 |
|
|
REAL, DIMENSION(klon) :: orogradT,xsi0 |
152 |
|
|
REAL, DIMENSION (klon,klev) :: zlay |
153 |
|
|
REAL :: Lvs, temp0 |
154 |
|
|
|
155 |
|
|
|
156 |
|
|
! Calculation of the near-surface temperature gradient along the topography |
157 |
|
|
!-------------------------------------------------------------------------- |
158 |
|
|
|
159 |
|
|
! at the moment, we fix it at a constant value (moist adiab. lapse rate) |
160 |
|
|
|
161 |
|
|
orogradT(:)=-6.5/1000. ! K/m |
162 |
|
|
|
163 |
|
|
! Calculation of near-surface surface ratqs |
164 |
|
|
!------------------------------------------- |
165 |
|
|
|
166 |
|
|
DO i=1,klon |
167 |
|
|
temp0=temp(i,1) |
168 |
|
|
IF (temp0 .LT. RTT) THEN |
169 |
|
|
Lvs=RLSTT |
170 |
|
|
ELSE |
171 |
|
|
Lvs=RLVTT |
172 |
|
|
ENDIF |
173 |
|
|
xsi0(i)=zstd(i)*ABS(orogradT(i))*Lvs/temp0/temp0/RV |
174 |
|
|
ratqs_oro(i,1)=xsi0(i) |
175 |
|
|
END DO |
176 |
|
|
|
177 |
|
|
! Vertical profile of ratqs assuming an exponential decrease with height |
178 |
|
|
!------------------------------------------------------------------------ |
179 |
|
|
|
180 |
|
|
! calculation of geop. height AGL |
181 |
|
|
zlay(:,1)= RD*temp(:,1)/(0.5*(paprs(:,1)+pplay(:,1))) & |
182 |
|
|
*(paprs(:,1)-pplay(:,1))/RG |
183 |
|
|
|
184 |
|
|
DO k=2,klev |
185 |
|
|
DO i = 1, klon |
186 |
|
|
zlay(i,k)= zlay(i,k-1)+RD*0.5*(temp(i,k-1)+temp(i,k)) & |
187 |
|
|
/paprs(i,k)*(pplay(i,k-1)-pplay(i,k))/RG |
188 |
|
|
|
189 |
|
|
ratqs_oro(i,k)=MAX(0.0,pctsrf(i,is_ter)*xsi0(i)*exp(-zlay(i,k)/MAX(zstd(i),1.))) |
190 |
|
|
END DO |
191 |
|
|
END DO |
192 |
|
|
|
193 |
|
|
|
194 |
|
|
|
195 |
|
|
|
196 |
|
|
END SUBROUTINE calcratqs_oro |
197 |
|
|
|
198 |
|
|
!============================================= |
199 |
|
|
|
200 |
|
|
SUBROUTINE calcratqs_hetero(klon,klev,t2m,q2m,temp,q,pplay,paprs,ratqs_hetero) |
201 |
|
|
|
202 |
|
|
! Etienne Vignon, November 2021 |
203 |
|
|
! Effect of subgrid surface heterogeneities on ratqs |
204 |
|
|
|
205 |
|
|
USE phys_local_var_mod, ONLY: s_pblh |
206 |
|
|
USE phys_state_var_mod, ONLY: pctsrf |
207 |
|
|
USE indice_sol_mod |
208 |
|
|
USE lscp_tools_mod, ONLY: CALC_QSAT_ECMWF |
209 |
|
|
|
210 |
|
|
IMPLICIT NONE |
211 |
|
|
|
212 |
|
|
include "YOMCST.h" |
213 |
|
|
|
214 |
|
|
! INPUTS |
215 |
|
|
|
216 |
|
|
|
217 |
|
|
INTEGER, INTENT(IN) :: klon ! number of horizontal grid points |
218 |
|
|
INTEGER, INTENT(IN) :: klev ! number of vertical layers |
219 |
|
|
REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: t2m ! 2m temperature for each tile [K] |
220 |
|
|
REAL, DIMENSION(klon,nbsrf), INTENT(IN) :: q2m ! 2m specific humidity for each tile [kg/kg] |
221 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! air temperature [K] |
222 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: q ! specific humidity [kg/kg] |
223 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! air pressure, layer's center [Pa] |
224 |
|
|
REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! air pressure, lower inteface [Pa] |
225 |
|
|
|
226 |
|
|
! OUTPUTS |
227 |
|
|
|
228 |
|
|
REAL, DIMENSION(klon,klev), INTENT(out) :: ratqs_hetero ! ratsq profile due to surface heterogeneities |
229 |
|
|
|
230 |
|
|
|
231 |
|
|
INTEGER :: i,k,nsrf |
232 |
|
|
REAL, DIMENSION(klon) :: xsi0, ratiom, qsat2m, dqsatdT |
233 |
|
|
REAL, DIMENSION (klon,klev) :: zlay |
234 |
|
|
|
235 |
|
|
|
236 |
|
|
|
237 |
|
|
! Calculation of near-surface surface ratqs |
238 |
|
|
!------------------------------------------- |
239 |
|
|
|
240 |
|
|
|
241 |
|
|
ratiom(:)=0. |
242 |
|
|
xsi0(:)=0. |
243 |
|
|
|
244 |
|
|
DO nsrf=1,nbsrf |
245 |
|
|
CALL CALC_QSAT_ECMWF(klon,t2m(:,nsrf),q2m(:,nsrf),paprs(:,1),RTT,0,.false.,qsat2m,dqsatdT) |
246 |
|
|
ratiom(:)=ratiom(:)+pctsrf(:,nsrf)*(q2m(:,nsrf)/qsat2m(:)) |
247 |
|
|
xsi0(:)=xsi0(:)+pctsrf(:,nsrf)*((q2m(:,nsrf)/qsat2m(:)-ratiom(:))**2) |
248 |
|
|
END DO |
249 |
|
|
|
250 |
|
|
xsi0(:)=sqrt(xsi0(:))/(ratiom(:)+1E-6) |
251 |
|
|
|
252 |
|
|
|
253 |
|
|
|
254 |
|
|
! Vertical profile of ratqs assuming an exponential decrease with height |
255 |
|
|
!------------------------------------------------------------------------ |
256 |
|
|
|
257 |
|
|
! calculation of geop. height AGL |
258 |
|
|
|
259 |
|
|
zlay(:,1)= RD*temp(:,1)/(0.5*(paprs(:,1)+pplay(:,1))) & |
260 |
|
|
*(paprs(:,1)-pplay(:,1))/RG |
261 |
|
|
ratqs_hetero(:,1)=xsi0(:) |
262 |
|
|
|
263 |
|
|
DO k=2,klev |
264 |
|
|
DO i = 1, klon |
265 |
|
|
zlay(i,k)= zlay(i,k-1)+RD*0.5*(temp(i,k-1)+temp(i,k)) & |
266 |
|
|
/paprs(i,k)*(pplay(i,k-1)-pplay(i,k))/RG |
267 |
|
|
|
268 |
|
|
ratqs_hetero(i,k)=MAX(xsi0(i)*exp(-zlay(i,k)/(s_pblh(i)+1.0)),0.0) |
269 |
|
|
END DO |
270 |
|
|
END DO |
271 |
|
|
|
272 |
|
|
END SUBROUTINE calcratqs_hetero |
273 |
|
|
|
274 |
|
|
!============================================= |
275 |
|
|
|
276 |
|
|
SUBROUTINE calcratqs_tke(klon,klev,pdtphys,temp,q,qsat,pplay,paprs,tke,tke_dissip,lmix,wprime,ratqs_tke) |
277 |
|
|
|
278 |
|
|
! References: |
279 |
|
|
! |
280 |
|
|
! Etienne Vignon: effect of subgrid turbulence on ratqs |
281 |
|
|
! |
282 |
|
|
! Field, P.R., Hill, A., Furtado, K., Korolev, A., 2014b. Mixed-phase clouds in a turbulent environment. Part |
283 |
|
|
! 2: analytic treatment. Q. J. R. Meteorol. Soc. 21, 2651–2663. https://doi.org/10.1002/qj.2175. |
284 |
|
|
! |
285 |
|
|
! Furtado, K., Field, P.R., Boutle, I.A., Morcrette, C.R., Wilkinson, J., 2016. A physically-based, subgrid |
286 |
|
|
! parametrization for the production and maintenance of mixed-phase clouds in a general circulation |
287 |
|
|
! model. J. Atmos. Sci. 73, 279–291. https://doi.org/10.1175/JAS-D-15-0021. |
288 |
|
|
|
289 |
|
|
USE phys_local_var_mod, ONLY: omega |
290 |
|
|
|
291 |
|
|
IMPLICIT NONE |
292 |
|
|
|
293 |
|
|
! INPUTS |
294 |
|
|
|
295 |
|
|
INTEGER, INTENT(IN) :: klon ! number of horizontal grid points |
296 |
|
|
INTEGER, INTENT(IN) :: klev ! number of vertical layers |
297 |
|
|
REAL, INTENT(IN) :: pdtphys ! physics time step [s] |
298 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: temp ! air temperature [K] |
299 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: q ! specific humidity [kg/kg] |
300 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: qsat ! saturation specific humidity [kg/kg] |
301 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! air pressure, layer's center [Pa] |
302 |
|
|
REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! air pressure, lower inteface [Pa] |
303 |
|
|
REAL, DIMENSION(klon,klev+1), INTENT(IN) :: tke ! Turbulent Kinetic Energy [m2/s2] |
304 |
|
|
REAL, DIMENSION(klon,klev+1), INTENT(IN) :: tke_dissip ! Turbulent Kinetic Energy Dissipation rate [m2/s3] |
305 |
|
|
REAL, DIMENSION(klon,klev+1), INTENT(IN) :: lmix ! Turbulent mixing length |
306 |
|
|
REAL, DIMENSION(klon,klev+1), INTENT(IN) :: wprime ! Turbulent vertical velocity scale [m/s] |
307 |
|
|
|
308 |
|
|
! OUTPUTS |
309 |
|
|
|
310 |
|
|
REAL, DIMENSION(klon,klev), INTENT(out) :: ratqs_tke ! ratsq profile due to subgrid TKE |
311 |
|
|
|
312 |
|
|
! LOCAL |
313 |
|
|
INTEGER :: i, k |
314 |
|
|
REAL :: AA, DD, NW, AAprime, VARLOG,rho,Lvs,taue,lhomo,dissmin,maxvarlog |
315 |
|
|
REAL, DIMENSION(klon,klev) :: sigmaw,w |
316 |
|
|
REAL, PARAMETER :: C0=10.0 |
317 |
|
|
REAL, PARAMETER :: lmin=0.001 |
318 |
|
|
REAL, PARAMETER :: ratqsmin=1E-6 |
319 |
|
|
REAL, PARAMETER :: ratqsmax=0.5 |
320 |
|
|
|
321 |
|
|
|
322 |
|
|
! Calculation of large scale and turbulent vertical velocities |
323 |
|
|
!--------------------------------------------------------------- |
324 |
|
|
|
325 |
|
|
DO k=1,klev |
326 |
|
|
DO i=1,klon |
327 |
|
|
rho=pplay(i,k)/temp(i,k)/RD |
328 |
|
|
w(i,k)=-rho*RG*omega(i,k) |
329 |
|
|
sigmaw(i,k)=0.5*(wprime(i,k+1)+wprime(i,k)) ! turbulent vertical velocity at the middle of model layers. |
330 |
|
|
END DO |
331 |
|
|
END DO |
332 |
|
|
|
333 |
|
|
! Calculation of ratqs |
334 |
|
|
!--------------------------------------------------------------- |
335 |
|
|
ratqs_tke(:,1)=ratqsmin ! set to a very low value to avoid division by 0 in order parts |
336 |
|
|
! of the code |
337 |
|
|
DO k=2,klev ! we start from second model level since TKE is not defined at k=1 |
338 |
|
|
DO i=1,klon |
339 |
|
|
|
340 |
|
|
IF (temp(i,k) .LT. RTT) THEN |
341 |
|
|
Lvs=RLSTT |
342 |
|
|
ELSE |
343 |
|
|
Lvs=RLVTT |
344 |
|
|
ENDIF |
345 |
|
|
dissmin=0.01*(0.5*(tke(i,k)+tke(i,k+1))/pdtphys) |
346 |
|
|
maxvarlog=LOG(1.0+ratqsmax**2)! to prevent ratqs from exceeding an arbitrary threshold value |
347 |
|
|
AA=RG*(Lvs/(RCPD*temp(i,k)*temp(i,k)*RV) - 1./(RD*temp(i,k))) |
348 |
|
|
lhomo=MAX(0.5*(lmix(i,k)+lmix(i,k+1)),lmin) |
349 |
|
|
taue=(lhomo*lhomo/MAX(0.5*(tke_dissip(i,k)+tke_dissip(i,k+1)),dissmin))**(1./3) ! Fields et al. 2014 |
350 |
|
|
DD=1.0/taue |
351 |
|
|
NW=(sigmaw(i,k)**2)*SQRT(2./(C0*MAX(0.5*(tke_dissip(i,k)+tke_dissip(i,k+1)),dissmin))) |
352 |
|
|
AAprime=AA*NW |
353 |
|
|
VARLOG=AAprime/2./DD |
354 |
|
|
VARLOG=MIN(VARLOG,maxvarlog) |
355 |
|
|
ratqs_tke(i,k)=SQRT(MAX(EXP(VARLOG)-1.0,ratqsmin)) |
356 |
|
|
END DO |
357 |
|
|
END DO |
358 |
|
|
END SUBROUTINE calcratqs_tke |
359 |
|
|
|
360 |
|
|
|
361 |
|
|
|
362 |
|
|
|
363 |
|
|
|
364 |
|
|
END MODULE calcratqs_multi_mod |