Newer
Older
!MNH_LIC Copyright 2013-2023 CNRS, Meteo-France and Universite Paul Sabatier
!MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
!MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
!MNH_LIC for details. version 1.
!-----------------------------------------------------------------
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
!
MODULE MODE_ELEC_BEARD_EFFECT
!
IMPLICIT NONE
CONTAINS
!
! ###################################################################
SUBROUTINE ELEC_BEARD_EFFECT(D, KID, OSEDIM, PT, PRHODREF, &
PRX, PQX, PEFIELDW, PLBDA, PBEARDCOEF)
! ####################################################################
!
!! PURPOSE
!! -------
!! The purpose of this routine is to compute the effect of the electric field
!! on the terminal velocity of hydrometeors.
!!
!! METHOD
!! ------
!! From Beard, K. V., 1980: The Effects of Altitude and Electrical Force on
!! the Terminal Velocity of Hydrometeors. J. Atmos. Sci., 37, 1363–1374,
!! https://doi.org/10.1175/1520-0469(1980)037<1363:TEOAAE>2.0.CO;2.
!!
!! AUTHOR
!! ------
!! J.-P. Pinty * LAERO *
!! C. Barthe * LAERO *
!!
!! MODIFICATIONS
!! -------------
!! Original 21/08/2013 first coded in rain_ice_elec
!! C. Barthe 01/06/2023 : externalize the code to use it with ICE3 and LIMA
!! C. Barthe 08/06/2023 : correction by 10-5 of the dynamic viscosity of air
!! (unecessary for eta0/eta but necessary for Re0)
!
!-------------------------------------------------------------------------------
!
!* 0. DECLARATIONS
! ------------
!
USE MODD_CST, ONLY: XG, XRD, XP00, XTT
USE MODD_DIMPHYEX, ONLY: DIMPHYEX_t
USE MODD_ELEC_DESCR, ONLY: XRTMIN_ELEC
USE MODD_PARAMETERS, ONLY: JPVEXT
USE MODD_PARAM_LIMA, ONLY: XALPHAC_L=>XALPHAC, XNUC_L=>XNUC, XALPHAR_L=>XALPHAR, XNUR_L=>XNUR, &
XALPHAI_L=>XALPHAI, XNUI_L=>XNUI, XALPHAS_L=>XALPHAS, XNUS_L=>XNUS, &
XALPHAG_L=>XALPHAG, XNUG_L=>XNUG, &
XCEXVT_L=>XCEXVT
USE MODD_PARAM_LIMA_COLD, ONLY: XBI_L=>XBI, XC_I_L=>XC_I, XDI_L=>XDI, &
XBS_L=>XBS, XCS_L=>XCS, XDS_L=>XDS
USE MODD_PARAM_LIMA_MIXED,ONLY: XBG_L=>XBG, XCG_L=>XCG, XDG_L=>XDG, &
XBH_L=>XBH, XCH_L=>XCH, XDH_L=>XDH, &
XALPHAH_L=>XALPHAH, XNUH_L=>XNUH
USE MODD_PARAM_LIMA_WARM, ONLY: XBR_L=>XBR, XCR_L=>XCR, XDR_L=>XDR, &
XBC_L=>XBC, XCC_L=>XCC, XDC_L=>XDC
USE MODD_PARAM_n, ONLY: CCLOUD
RODIER Quentin
committed
USE MODD_RAIN_ICE_DESCR_n,ONLY: XALPHAC_I=>XALPHAC, XNUC_I=>XNUC, XALPHAR_I=>XALPHAR, XNUR_I=>XNUR, &
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
XALPHAI_I=>XALPHAI, XNUI_I=>XNUI, XALPHAS_I=>XALPHAS, XNUS_I=>XNUS, &
XALPHAG_I=>XALPHAG, XNUG_I=>XNUG, XALPHAH_I=>XALPHAH, XNUH_I=>XNUH, &
XBC_I=>XBC, XCC_I=>XCC, XDC_I=>XDC, &
XBR_I=>XBR, XCR_I=>XCR, XDR_I=>XDR, &
XBI_I=>XBI, XC_I_I=>XC_I, XDI_I=>XDI, &
XBS_I=>XBS, XCS_I=>XCS, XDS_I=>XDS, &
XBG_I=>XBG, XCG_I=>XCG, XDG_I=>XDG, &
XBH_I=>XBH, XCH_I=>XCH, XDH_I=>XDH, &
XCEXVT_I=>XCEXVT
USE MODD_REF, ONLY: XTHVREFZ
!
USE MODI_MOMG
!
USE MODE_TOOLS_ll, ONLY: GET_INDICE_ll
!
IMPLICIT NONE
!
!* 0.1 Declarations of dummy arguments
!
TYPE(DIMPHYEX_t), INTENT(IN) :: D
INTEGER, INTENT(IN) :: KID ! Hydrometeor ID
LOGICAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: OSEDIM ! if T, compute the sedim. proc.
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRHODREF ! Reference density
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PT ! Temperature
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PRX ! m.r. source
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PQX ! Elec. charge density source
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PEFIELDW ! Vertical component of the electric field
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(IN) :: PLBDA ! Slope param. of the distribution
REAL, DIMENSION(D%NIJT,D%NKT), INTENT(INOUT) :: PBEARDCOEF ! Beard coefficient
!
!* 0.2 Declarations of local variables
!
INTEGER :: JIJ, JK ! loop indexes
INTEGER :: IIJB, IIJE, IKTB, IKTE
REAL :: ZCEXVT, ZBX, ZCX, ZDX, ZALPHAX, ZNUX
REAL :: ZRE0
REAL :: ZETA0
REAL :: ZVX
REAL :: ZK
REAL :: ZCOR00, ZRHO00
REAL :: ZT ! Temperature (C)
REAL :: ZCOR ! To remove the Foote-duToit correction
REAL :: ZF0, ZF1 ! Coef. in Beard's equation
real, dimension(D%NIJT,D%NKT) :: zreynolds
!
!-------------------------------------------------------------------------------
!
!* 1. COMPUTE USEFULL PARAMETERS
! --------------------------
!
IKTB = D%NKTB
IKTE = D%NKTE
IIJB = D%NIJB
IIJE = D%NIJE
!
!* 1.1 Select the right parameters
! --> depend on the microphysics scheme and the hydrometeor species
!
IF (CCLOUD(1:3) == 'ICE') THEN
ZCEXVT = XCEXVT_I
!
IF (KID == 2) THEN
ZBX = XBC_I
ZCX = XCC_I
ZDX = XDC_I
ZALPHAX = XALPHAC_I
ZNUX = XNUC_I
ELSE IF (KID == 3) THEN
ZBX = XBR_I
ZCX = XCR_I
ZDX = XDR_I
ZALPHAX = XALPHAR_I
ZNUX = XNUR_I
ELSE IF (KID == 4) THEN
! values for columns are used to be consistent with the McF&H formula
ZBX = 1.7
ZCX = 2.1E5
ZDX = 1.585
ZALPHAX = XALPHAI_I
ZNUX = XNUI_I
ELSE IF (KID == 5) THEN
ZBX = XBS_I
ZCX = XCS_I
ZDX = XDS_I
ZALPHAX = XALPHAS_I
ZNUX = XNUS_I
ELSE IF (KID == 6) THEN
ZBX = XBG_I
ZCX = XCG_I
ZDX = XDG_I
ZALPHAX = XALPHAG_I
ZNUX = XNUG_I
ELSE IF (KID == 7) THEN
ZBX = XBH_I
ZCX = XCH_I
ZDX = XDH_I
ZALPHAX = XALPHAH_I
ZNUX = XNUH_I
END IF
ELSE IF (CCLOUD == 'LIMA') THEN
ZCEXVT = XCEXVT_L
!
IF (KID == 2) THEN
ZBX = XBC_L
ZCX = XCC_L
ZDX = XDC_L
ZALPHAX = XALPHAC_L
ZNUX = XNUC_L
ELSE IF (KID == 3) THEN
ZBX = XBR_L
ZCX = XCR_L
ZDX = XDR_L
ZALPHAX = XALPHAR_L
ZNUX = XNUR_L
ELSE IF (KID == 4) THEN
ZBX = 1.7
ZCX = 2.1E5
ZDX = 1.585
ZALPHAX = XALPHAI_L
ZNUX = XNUI_L
ELSE IF (KID == 5) THEN
ZBX = XBS_L
ZCX = XCS_L
ZDX = XDS_L
ZALPHAX = XALPHAS_L
ZNUX = XNUS_L
ELSE IF (KID == 6) THEN
ZBX = XBG_L
ZCX = XCG_L
ZDX = XDG_L
ZALPHAX = XALPHAG_L
ZNUX = XNUG_L
ELSE IF (KID == 7) THEN
ZBX = XBH_L
ZCX = XCH_L
ZDX = XDH_L
ZALPHAX = XALPHAH_L
ZNUX = XNUH_L
END IF
!
END IF
!
!* 1.2 Parameters from Table 1 in Beard (1980)
!
! Reference value of the dynamic viscosity of air
ZETA0 = (1.718E-5 + 0.0049E-5 * (XTHVREFZ(IKTB) - XTT))
!
ZRHO00 = XP00 / (XRD * XTHVREFZ(IKTB))
ZCOR00 = ZRHO00**ZCEXVT
!
! (rho_0 / eta_0) * (v * lambda^d)
ZVX = (ZRHO00 / ZETA0) * ZCX * MOMG(ZALPHAX,ZNUX,ZBX+ZDX) / MOMG(ZALPHAX,ZNUX,ZBX)
!
!-------------------------------------------------------------------------------
!
!* 2. COMPUTE THE VELOCITY ADJUSTMENT FACTOR
! --------------------------------------
!
zreynolds(:,:) = 0.
PBEARDCOEF(:,:) = 1.0
!
DO JK = IKTB, IKTE
DO JIJ = IIJB, IIJE
!++cb++ 09/06/23 on n'applique l'effet Beard que pour les points ou le rapport de melange est
! suffisamment eleve pour eviter que qE >> mg => coef de Beard tres eleve !
! Ce pb intervient avec ICE3 pour lequel xrtmin est tres bas par rapport a LIMA.
IF (OSEDIM(JIJ,JK) .AND. PRX(JIJ,JK) .GT. XRTMIN_ELEC(KID) .AND. PLBDA(JIJ,JK) .GT. 0.) THEN
!--cb--
! Temperature K --> C
ZT = PT(JIJ,JK) - XTT
!
! Pre-factor of f_0
IF (ZT >= 0.0) THEN
ZF0 = ZETA0 / (1.718E-5 + 0.0049E-5 * ZT)
ELSE
ZF0 = ZETA0 / (1.718E-5 + 0.0049E-5 * ZT - 1.2E-10 * ZT * ZT)
END IF
!
! Pre-factor of f_infty
ZF1 = SQRT(ZRHO00/PRHODREF(JIJ,JK))
!
! compute (1 - K) = 1 - qE/mg
ZK = 1. - PQX(JIJ,JK) * PEFIELDW(JIJ,JK) / (PRX(JIJ,JK) * XG)
!
! Hyp : K_0 ~ 0
! Hyp : si qE > mg, K > 1
IF (ZK <= 0.0) THEN
PBEARDCOEF(JIJ,JK) = 0. ! levitation
ELSE
! Reynolds number
ZRE0 = ZVX / PLBDA(JIJ,JK)**(1.+ZDX)
zreynolds(jij,jk) = zre0
IF (ZRE0 <= 0.2) THEN
PBEARDCOEF(JIJ,JK) = ZF0 * ZK
ELSE IF (ZRE0 >= 1000.) THEN
PBEARDCOEF(JIJ,JK) = ZF1 * SQRT(ZK)
ELSE
PBEARDCOEF(JIJ,JK) = ZF0 * ZK + &
(ZF1 * SQRT(ZK) - ZF0 * ZK) * &
(1.61 + LOG(ZRE0)) / 8.52
END IF
! remove the Foote-duToit correction
ZCOR = (PRHODREF(JIJ,JK) / ZRHO00)**ZCEXVT
PBEARDCOEF(JIJ,JK) = PBEARDCOEF(JIJ,JK) * ZCOR
END IF
ELSE
PBEARDCOEF(JIJ,JK) = 1.0 ! No "Beard" effect
END IF
END DO
END DO
!
!-------------------------------------------------------------------------------
!
END SUBROUTINE ELEC_BEARD_EFFECT
END MODULE MODE_ELEC_BEARD_EFFECT