1 |
|
|
MODULE lscp_mod |
2 |
|
|
|
3 |
|
|
IMPLICIT NONE |
4 |
|
|
|
5 |
|
|
CONTAINS |
6 |
|
|
|
7 |
|
|
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
8 |
|
288 |
SUBROUTINE lscp(klon,klev,dtime,missing_val, & |
9 |
|
288 |
paprs,pplay,t,q,ptconv,ratqs, & |
10 |
|
288 |
d_t, d_q, d_ql, d_qi, rneb, rneblsvol, rneb_seri, & |
11 |
|
|
pfraclr,pfracld, & |
12 |
|
288 |
radocond, radicefrac, rain, snow, & |
13 |
|
|
frac_impa, frac_nucl, beta, & |
14 |
|
288 |
prfl, psfl, rhcl, zqta, fraca, & |
15 |
|
|
ztv, zpspsk, ztla, zthl, iflag_cld_th, & |
16 |
|
|
iflag_ice_thermo, ok_ice_sursat, qsatl, qsats, & |
17 |
|
|
distcltop, & |
18 |
|
|
qclr, qcld, qss, qvc, rnebclr, rnebss, gamma_ss, & |
19 |
|
|
Tcontr, qcontr, qcontr2, fcontrN, fcontrP) |
20 |
|
|
|
21 |
|
|
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
22 |
|
|
! Authors: Z.X. Li (LMD), J-L Dufresne (LMD), C. Rio (LMD), J-Y Grandpeix (LMD) |
23 |
|
|
! A. JAM (LMD), J-B Madeleine (LMD), E. Vignon (LMD), L. Touzze-Peiffert (LMD) |
24 |
|
|
!-------------------------------------------------------------------------------- |
25 |
|
|
! Date: 01/2021 |
26 |
|
|
!-------------------------------------------------------------------------------- |
27 |
|
|
! Aim: Large Scale Clouds and Precipitation (LSCP) |
28 |
|
|
! |
29 |
|
|
! This code is a new version of the fisrtilp.F90 routine, which is itself a |
30 |
|
|
! merge of 'first' (superrsaturation physics, P. LeVan K. Laval) |
31 |
|
|
! and 'ilp' (il pleut, L. Li) |
32 |
|
|
! |
33 |
|
|
! Compared to the original fisrtilp code, lscp |
34 |
|
|
! -> assumes thermcep = .TRUE. all the time (fisrtilp inconsistent when .FALSE.) |
35 |
|
|
! -> consider always precipitation thermalisation (fl_cor_ebil>0) |
36 |
|
|
! -> option iflag_fisrtilp_qsat<0 no longer possible (qsat does not evolve with T) |
37 |
|
|
! -> option "oldbug" by JYG has been removed |
38 |
|
|
! -> iflag_t_glace >0 always |
39 |
|
|
! -> the 'all or nothing' cloud approach is no longer available (cpartiel=T always) |
40 |
|
|
! -> rectangular distribution from L. Li no longer available |
41 |
|
|
! -> We always account for the Wegener-Findeisen-Bergeron process (iflag_bergeron = 2 in fisrt) |
42 |
|
|
!-------------------------------------------------------------------------------- |
43 |
|
|
! References: |
44 |
|
|
! |
45 |
|
|
! - Bony, S., & Emanuel, K. A. 2001, JAS, doi: 10.1175/1520-0469(2001)058<3158:APOTCA>2.0.CO;2 |
46 |
|
|
! - Hourdin et al. 2013, Clim Dyn, doi:10.1007/s00382-012-1343-y |
47 |
|
|
! - Jam et al. 2013, Boundary-Layer Meteorol, doi:10.1007/s10546-012-9789-3 |
48 |
|
|
! - Jouhaud, et al. 2018. JAMES, doi:10.1029/2018MS001379 |
49 |
|
|
! - Madeleine et al. 2020, JAMES, doi:10.1029/2020MS002046 |
50 |
|
|
! - Touzze-Peifert Ludo, PhD thesis, p117-124 |
51 |
|
|
! ------------------------------------------------------------------------------- |
52 |
|
|
! Code structure: |
53 |
|
|
! |
54 |
|
|
! P0> Thermalization of the precipitation coming from the overlying layer |
55 |
|
|
! P1> Evaporation of the precipitation (falling from the k+1 level) |
56 |
|
|
! P2> Cloud formation (at the k level) |
57 |
|
|
! P2.A.1> With the PDFs, calculation of cloud properties using the inital |
58 |
|
|
! values of T and Q |
59 |
|
|
! P2.A.2> Coupling between condensed water and temperature |
60 |
|
|
! P2.A.3> Calculation of final quantities associated with cloud formation |
61 |
|
|
! P2.B> Release of Latent heat after cloud formation |
62 |
|
|
! P3> Autoconversion to precipitation (k-level) |
63 |
|
|
! P4> Wet scavenging |
64 |
|
|
!------------------------------------------------------------------------------ |
65 |
|
|
! Some preliminary comments (JBM) : |
66 |
|
|
! |
67 |
|
|
! The cloud water that the radiation scheme sees is not the same that the cloud |
68 |
|
|
! water used in the physics and the dynamics |
69 |
|
|
! |
70 |
|
|
! During the autoconversion to precipitation (P3 step), radocond (cloud water used |
71 |
|
|
! by the radiation scheme) is calculated as an average of the water that remains |
72 |
|
|
! in the cloud during the precipitation and not the water remaining at the end |
73 |
|
|
! of the time step. The latter is used in the rest of the physics and advected |
74 |
|
|
! by the dynamics. |
75 |
|
|
! |
76 |
|
|
! In summary: |
77 |
|
|
! |
78 |
|
|
! Radiation: |
79 |
|
|
! xflwc(newmicro)+xfiwc(newmicro) = |
80 |
|
|
! radocond=lwcon(nc)+iwcon(nc) |
81 |
|
|
! |
82 |
|
|
! Notetheless, be aware of: |
83 |
|
|
! |
84 |
|
|
! radocond .NE. ocond(nc) |
85 |
|
|
! i.e.: |
86 |
|
|
! lwcon(nc)+iwcon(nc) .NE. ocond(nc) |
87 |
|
|
! but oliq+(ocond-oliq) .EQ. ocond |
88 |
|
|
! (which is not trivial) |
89 |
|
|
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
90 |
|
|
! |
91 |
|
|
USE print_control_mod, ONLY: prt_level, lunout |
92 |
|
|
USE cloudth_mod, ONLY : cloudth, cloudth_v3, cloudth_v6, cloudth_mpc |
93 |
|
|
USE lscp_tools_mod, ONLY : calc_qsat_ecmwf, icefrac_lscp, calc_gammasat |
94 |
|
|
USE lscp_tools_mod, ONLY : fallice_velocity, distance_to_cloud_top |
95 |
|
|
USE ice_sursat_mod, ONLY : ice_sursat |
96 |
|
|
|
97 |
|
|
USE lscp_ini_mod, ONLY : seuil_neb, niter_lscp, iflag_evap_prec, t_coup, DDT0, ztfondue, rain_int_min |
98 |
|
|
USE lscp_ini_mod, ONLY : iflag_mpc_bl, ok_radocond_snow, a_tr_sca, cld_expo_con, cld_expo_lsc |
99 |
|
|
USE lscp_ini_mod, ONLY : iflag_cloudth_vert, iflag_rain_incloud_vol, iflag_t_glace, t_glace_min |
100 |
|
|
USE lscp_ini_mod, ONLY : coef_eva, coef_eva_i,cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con |
101 |
|
|
USE lscp_ini_mod, ONLY : iflag_bergeron, iflag_fisrtilp_qsat, iflag_vice, cice_velo, dice_velo |
102 |
|
|
USE lscp_ini_mod, ONLY : iflag_autoconversion, ffallv_con, ffallv_lsc |
103 |
|
|
USE lscp_ini_mod, ONLY : RCPD, RLSTT, RLVTT, RLMLT, RVTMP2, RTT, RD, RG |
104 |
|
|
|
105 |
|
|
IMPLICIT NONE |
106 |
|
|
|
107 |
|
|
!=============================================================================== |
108 |
|
|
! VARIABLES DECLARATION |
109 |
|
|
!=============================================================================== |
110 |
|
|
|
111 |
|
|
! INPUT VARIABLES: |
112 |
|
|
!----------------- |
113 |
|
|
|
114 |
|
|
INTEGER, INTENT(IN) :: klon,klev ! number of horizontal grid points and vertical levels |
115 |
|
|
REAL, INTENT(IN) :: dtime ! time step [s] |
116 |
|
|
REAL, INTENT(IN) :: missing_val ! missing value for output |
117 |
|
|
|
118 |
|
|
REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! inter-layer pressure [Pa] |
119 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! mid-layer pressure [Pa] |
120 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: t ! temperature (K) |
121 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: q ! total specific humidity (= vapor phase) [kg/kg] |
122 |
|
|
INTEGER, INTENT(IN) :: iflag_cld_th ! flag that determines the distribution of convective clouds |
123 |
|
|
INTEGER, INTENT(IN) :: iflag_ice_thermo! flag to activate the ice thermodynamics |
124 |
|
|
! CR: if iflag_ice_thermo=2, only convection is active |
125 |
|
|
LOGICAL, INTENT(IN) :: ok_ice_sursat ! flag to determine if ice sursaturation is activated |
126 |
|
|
|
127 |
|
|
LOGICAL, DIMENSION(klon,klev), INTENT(IN) :: ptconv ! grid points where deep convection scheme is active |
128 |
|
|
|
129 |
|
|
!Inputs associated with thermal plumes |
130 |
|
|
|
131 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: ztv ! virtual potential temperature [K] |
132 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: zqta ! specific humidity within thermals [kg/kg] |
133 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: fraca ! fraction of thermals within the mesh [-] |
134 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: zpspsk ! exner potential (p/100000)**(R/cp) |
135 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: ztla ! liquid temperature within thermals [K] |
136 |
|
|
|
137 |
|
|
! INPUT/OUTPUT variables |
138 |
|
|
!------------------------ |
139 |
|
|
|
140 |
|
|
REAL, DIMENSION(klon,klev), INTENT(INOUT) :: zthl ! liquid potential temperature [K] |
141 |
|
|
REAL, DIMENSION(klon,klev), INTENT(INOUT):: ratqs ! function of pressure that sets the large-scale |
142 |
|
|
REAL, DIMENSION(klon,klev), INTENT(INOUT):: beta ! conversion rate of condensed water |
143 |
|
|
|
144 |
|
|
|
145 |
|
|
! Input sursaturation en glace |
146 |
|
|
REAL, DIMENSION(klon,klev), INTENT(INOUT):: rneb_seri ! fraction nuageuse en memoire |
147 |
|
|
|
148 |
|
|
! OUTPUT variables |
149 |
|
|
!----------------- |
150 |
|
|
|
151 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t ! temperature increment [K] |
152 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q ! specific humidity increment [kg/kg] |
153 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_ql ! liquid water increment [kg/kg] |
154 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_qi ! cloud ice mass increment [kg/kg] |
155 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: rneb ! cloud fraction [-] |
156 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: rneblsvol ! cloud fraction per unit volume [-] |
157 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfraclr ! precip fraction clear-sky part [-] |
158 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: pfracld ! precip fraction cloudy part [-] |
159 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: radocond ! condensed water used in the radiation scheme [kg/kg] |
160 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: radicefrac ! ice fraction of condensed water for radiation scheme |
161 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: rhcl ! clear-sky relative humidity [-] |
162 |
|
|
REAL, DIMENSION(klon), INTENT(OUT) :: rain ! surface large-scale rainfall [kg/s/m2] |
163 |
|
|
REAL, DIMENSION(klon), INTENT(OUT) :: snow ! surface large-scale snowfall [kg/s/m2] |
164 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsatl ! saturation specific humidity wrt liquid [kg/kg] |
165 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: qsats ! saturation specific humidity wrt ice [kg/kg] |
166 |
|
|
REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: prfl ! large-scale rainfall flux in the column [kg/s/m2] |
167 |
|
|
REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: psfl ! large-scale snowfall flux in the column [kg/s/m2] |
168 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: distcltop ! distance to cloud top [m] |
169 |
|
|
|
170 |
|
|
! fraction of aerosol scavenging through impaction and nucleation (for on-line) |
171 |
|
|
|
172 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_impa ! scavenging fraction due tu impaction [-] |
173 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_nucl ! scavenging fraction due tu nucleation [-] |
174 |
|
|
|
175 |
|
|
! for supersaturation and contrails parameterisation |
176 |
|
|
|
177 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: qclr ! specific total water content in clear sky region [kg/kg] |
178 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcld ! specific total water content in cloudy region [kg/kg] |
179 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: qss ! specific total water content in supersat region [kg/kg] |
180 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: qvc ! specific vapor content in clouds [kg/kg] |
181 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: rnebclr ! mesh fraction of clear sky [-] |
182 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: rnebss ! mesh fraction of ISSR [-] |
183 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: gamma_ss ! coefficient governing the ice nucleation RHi threshold [-] |
184 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: Tcontr ! threshold temperature for contrail formation [K] |
185 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcontr ! threshold humidity for contrail formation [kg/kg] |
186 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: qcontr2 ! // (2nd expression more consistent with LMDZ expression of q) |
187 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: fcontrN ! fraction of grid favourable to non-persistent contrails |
188 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: fcontrP ! fraction of grid favourable to persistent contrails |
189 |
|
|
|
190 |
|
|
|
191 |
|
|
! LOCAL VARIABLES: |
192 |
|
|
!---------------- |
193 |
|
|
|
194 |
|
576 |
REAL qsl(klon), qsi(klon) |
195 |
|
|
REAL zct, zcl,zexpo |
196 |
|
576 |
REAL ctot(klon,klev) |
197 |
|
576 |
REAL ctot_vol(klon,klev) |
198 |
|
576 |
REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 |
199 |
|
576 |
REAL zdqsdT_raw(klon) |
200 |
|
576 |
REAL gammasat(klon),dgammasatdt(klon) ! coefficient to make cold condensation at the correct RH and derivative wrt T |
201 |
|
576 |
REAL Tbef(klon),qlbef(klon),DT(klon),num,denom |
202 |
|
|
REAL cste |
203 |
|
576 |
REAL zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon) |
204 |
|
576 |
REAL Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon) |
205 |
|
|
REAL erf |
206 |
|
576 |
REAL qcloud(klon), icefrac_mpc(klon,klev), qincloud_mpc(klon) |
207 |
|
576 |
REAL zrfl(klon), zrfln(klon), zqev, zqevt |
208 |
|
576 |
REAL zifl(klon), zifln(klon), ziflprev(klon),zqev0,zqevi, zqevti |
209 |
|
576 |
REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon) |
210 |
|
576 |
REAL zoliql(klon), zoliqi(klon) |
211 |
|
576 |
REAL zt(klon) |
212 |
|
576 |
REAL zdz(klon),zrho(klon,klev),iwc(klon) |
213 |
|
576 |
REAL zchau,zfroi,zfice(klon),zneb(klon),znebprecip(klon) |
214 |
|
|
REAL zmelt,zrain,zsnow,zprecip |
215 |
|
576 |
REAL dzfice(klon) |
216 |
|
|
REAL zsolid |
217 |
|
576 |
REAL qtot(klon), qzero(klon) |
218 |
|
576 |
REAL dqsl(klon),dqsi(klon) |
219 |
|
|
REAL smallestreal |
220 |
|
|
! Variables for Bergeron process |
221 |
|
|
REAL zcp, coef1, DeltaT, Deltaq, Deltaqprecl |
222 |
|
576 |
REAL zqpreci(klon), zqprecl(klon) |
223 |
|
|
! Variables precipitation energy conservation |
224 |
|
576 |
REAL zmqc(klon) |
225 |
|
|
REAL zalpha_tr |
226 |
|
|
REAL zfrac_lessi |
227 |
|
576 |
REAL zprec_cond(klon) |
228 |
|
576 |
REAL zmair(klon), zcpair, zcpeau |
229 |
|
|
REAL zlh_solid(klon), zm_solid ! for liquid -> solid conversion |
230 |
|
576 |
REAL zrflclr(klon), zrflcld(klon) |
231 |
|
576 |
REAL d_zrfl_clr_cld(klon), d_zifl_clr_cld(klon) |
232 |
|
576 |
REAL d_zrfl_cld_clr(klon), d_zifl_cld_clr(klon) |
233 |
|
576 |
REAL ziflclr(klon), ziflcld(klon) |
234 |
|
576 |
REAL znebprecipclr(klon), znebprecipcld(klon) |
235 |
|
576 |
REAL tot_zneb(klon), tot_znebn(klon), d_tot_zneb(klon) |
236 |
|
576 |
REAL d_znebprecip_clr_cld(klon), d_znebprecip_cld_clr(klon) |
237 |
|
576 |
REAL velo(klon,klev), vr, ffallv |
238 |
|
|
REAL qlmpc, qimpc, rnebmpc |
239 |
|
576 |
REAL radocondi(klon,klev), radocondl(klon,klev) |
240 |
|
|
REAL effective_zneb |
241 |
|
576 |
REAL distcltop1D(klon) |
242 |
|
|
|
243 |
|
|
|
244 |
|
|
INTEGER i, k, n, kk |
245 |
|
576 |
INTEGER n_i(klon), iter |
246 |
|
|
INTEGER ncoreczq |
247 |
|
576 |
INTEGER mpc_bl_points(klon,klev) |
248 |
|
|
|
249 |
|
576 |
LOGICAL lognormale(klon) |
250 |
|
576 |
LOGICAL keepgoing(klon) |
251 |
|
|
LOGICAL,SAVE :: appel1er |
252 |
|
|
!$OMP THREADPRIVATE(appel1er) |
253 |
|
|
|
254 |
|
|
CHARACTER (len = 20) :: modname = 'lscp' |
255 |
|
|
CHARACTER (len = 80) :: abort_message |
256 |
|
|
|
257 |
|
|
DATA appel1er /.TRUE./ |
258 |
|
|
|
259 |
|
|
!=============================================================================== |
260 |
|
|
! INITIALISATION |
261 |
|
|
!=============================================================================== |
262 |
|
|
|
263 |
|
|
! Few initial checks |
264 |
|
|
|
265 |
✗✓ |
288 |
IF (iflag_fisrtilp_qsat .LT. 0) THEN |
266 |
|
|
abort_message = 'lscp cannot be used with iflag_fisrtilp<0' |
267 |
|
|
CALL abort_physic(modname,abort_message,1) |
268 |
|
|
ENDIF |
269 |
|
|
|
270 |
|
|
! Few initialisations |
271 |
|
|
|
272 |
✓✓ |
286560 |
znebprecip(:)=0.0 |
273 |
✓✓✓✓
|
11176128 |
ctot_vol(1:klon,1:klev)=0.0 |
274 |
✓✓✓✓
|
11176128 |
rneblsvol(1:klon,1:klev)=0.0 |
275 |
|
|
smallestreal=1.e-9 |
276 |
✓✓ |
286560 |
znebprecipclr(:)=0.0 |
277 |
✓✓ |
286560 |
znebprecipcld(:)=0.0 |
278 |
✓✓✓✓
|
11176128 |
mpc_bl_points(:,:)=0 |
279 |
|
|
|
280 |
✗✓ |
288 |
IF (prt_level>9) WRITE(lunout,*) 'NUAGES4 A. JAM' |
281 |
|
|
|
282 |
|
|
! AA for 'safety' reasons |
283 |
|
|
zalpha_tr = 0. |
284 |
|
|
zfrac_lessi = 0. |
285 |
✓✓✓✓
|
11176128 |
beta(:,:)= 0. |
286 |
|
|
|
287 |
|
|
! Initialisation of variables: |
288 |
|
|
|
289 |
✓✓✓✓
|
11462688 |
prfl(:,:) = 0.0 |
290 |
✓✓✓✓
|
11462688 |
psfl(:,:) = 0.0 |
291 |
✓✓✓✓
|
11176128 |
d_t(:,:) = 0.0 |
292 |
✓✓✓✓
|
11176128 |
d_q(:,:) = 0.0 |
293 |
✓✓✓✓
|
11176128 |
d_ql(:,:) = 0.0 |
294 |
✓✓✓✓
|
11176128 |
d_qi(:,:) = 0.0 |
295 |
✓✓✓✓
|
11176128 |
rneb(:,:) = 0.0 |
296 |
✓✓✓✓
|
11176128 |
pfraclr(:,:)=0.0 |
297 |
✓✓✓✓
|
11176128 |
pfracld(:,:)=0.0 |
298 |
✓✓✓✓
|
11176128 |
radocond(:,:) = 0.0 |
299 |
✓✓✓✓
|
11176128 |
radicefrac(:,:) = 0.0 |
300 |
✓✓✓✓
|
11176128 |
frac_nucl(:,:) = 1.0 |
301 |
✓✓✓✓
|
11176128 |
frac_impa(:,:) = 1.0 |
302 |
✓✓ |
286560 |
rain(:) = 0.0 |
303 |
✓✓ |
286560 |
snow(:) = 0.0 |
304 |
✓✓ |
286560 |
zoliq(:)=0.0 |
305 |
✓✓ |
286560 |
zfice(:)=0.0 |
306 |
✓✓ |
286560 |
dzfice(:)=0.0 |
307 |
✓✓ |
286560 |
zqprecl(:)=0.0 |
308 |
✓✓ |
286560 |
zqpreci(:)=0.0 |
309 |
✓✓ |
286560 |
zrfl(:) = 0.0 |
310 |
✓✓ |
286560 |
zifl(:) = 0.0 |
311 |
✓✓ |
286560 |
ziflprev(:)=0.0 |
312 |
✓✓ |
286560 |
zneb(:) = seuil_neb |
313 |
✓✓ |
286560 |
zrflclr(:) = 0.0 |
314 |
✓✓ |
286560 |
ziflclr(:) = 0.0 |
315 |
✓✓ |
286560 |
zrflcld(:) = 0.0 |
316 |
✓✓ |
286560 |
ziflcld(:) = 0.0 |
317 |
✓✓ |
286560 |
tot_zneb(:) = 0.0 |
318 |
✓✓ |
286560 |
tot_znebn(:) = 0.0 |
319 |
✓✓ |
286560 |
d_tot_zneb(:) = 0.0 |
320 |
✓✓ |
286560 |
qzero(:) = 0.0 |
321 |
✓✓ |
286560 |
distcltop1D(:)=0.0 |
322 |
|
|
|
323 |
|
|
!--ice supersaturation |
324 |
✓✓✓✓
|
11176128 |
gamma_ss(:,:) = 1.0 |
325 |
✓✓✓✓
|
11176128 |
qss(:,:) = 0.0 |
326 |
✓✓✓✓
|
11176128 |
rnebss(:,:) = 0.0 |
327 |
✓✓✓✓
|
11176128 |
Tcontr(:,:) = missing_val |
328 |
✓✓✓✓
|
11176128 |
qcontr(:,:) = missing_val |
329 |
✓✓✓✓
|
11176128 |
qcontr2(:,:) = missing_val |
330 |
✓✓✓✓
|
11176128 |
fcontrN(:,:) = 0.0 |
331 |
✓✓✓✓
|
11176128 |
fcontrP(:,:) = 0.0 |
332 |
|
|
|
333 |
|
|
!c_iso: variable initialisation for iso |
334 |
|
|
|
335 |
|
|
|
336 |
|
|
!=============================================================================== |
337 |
|
|
! BEGINNING OF VERTICAL LOOP FROM TOP TO BOTTOM |
338 |
|
|
!=============================================================================== |
339 |
|
|
|
340 |
|
288 |
ncoreczq=0 |
341 |
|
|
|
342 |
✓✓ |
11520 |
DO k = klev, 1, -1 |
343 |
|
|
|
344 |
|
|
! Initialisation temperature and specific humidity |
345 |
✓✓ |
11175840 |
DO i = 1, klon |
346 |
|
11164608 |
zt(i)=t(i,k) |
347 |
|
11175840 |
zq(i)=q(i,k) |
348 |
|
|
!c_iso init of iso |
349 |
|
|
ENDDO |
350 |
|
|
|
351 |
|
|
! -------------------------------------------------------------------- |
352 |
|
|
! P0> Thermalization of precipitation falling from the overlying layer |
353 |
|
|
! -------------------------------------------------------------------- |
354 |
|
|
! Computes air temperature variation due to enthalpy transported by |
355 |
|
|
! precipitation. Precipitation is then thermalized with the air in the |
356 |
|
|
! layer. |
357 |
|
|
! The precipitation should remain thermalized throughout the different |
358 |
|
|
! thermodynamical transformations. |
359 |
|
|
! The corresponding water mass should |
360 |
|
|
! be added when calculating the layer's enthalpy change with |
361 |
|
|
! temperature |
362 |
|
|
! See lmdzpedia page todoan |
363 |
|
|
! todoan: check consistency with ice phase |
364 |
|
|
! todoan: understand why several steps |
365 |
|
|
! --------------------------------------------------------------------- |
366 |
|
|
|
367 |
✓✓ |
11232 |
IF (k.LE.klev-1) THEN |
368 |
|
|
|
369 |
✓✓ |
10889280 |
DO i = 1, klon |
370 |
|
|
|
371 |
|
10878336 |
zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG |
372 |
|
|
! no condensed water so cp=cp(vapor+dry air) |
373 |
|
|
! RVTMP2=rcpv/rcpd-1 |
374 |
|
10878336 |
zcpair=RCPD*(1.0+RVTMP2*zq(i)) |
375 |
|
10878336 |
zcpeau=RCPD*RVTMP2 |
376 |
|
|
|
377 |
|
|
! zmqc: precipitation mass that has to be thermalized with |
378 |
|
|
! layer's air so that precipitation at the ground has the |
379 |
|
|
! same temperature as the lowermost layer |
380 |
|
10878336 |
zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair(i) |
381 |
|
|
! t(i,k+1)+d_t(i,k+1): new temperature of the overlying layer |
382 |
|
|
zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) & |
383 |
|
10889280 |
/ (zcpair + zmqc(i)*zcpeau) |
384 |
|
|
|
385 |
|
|
ENDDO |
386 |
|
|
|
387 |
|
|
ELSE |
388 |
|
|
|
389 |
✓✓ |
286560 |
DO i = 1, klon |
390 |
|
286272 |
zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG |
391 |
|
286560 |
zmqc(i) = 0. |
392 |
|
|
ENDDO |
393 |
|
|
|
394 |
|
|
ENDIF |
395 |
|
|
|
396 |
|
|
! -------------------------------------------------------------------- |
397 |
|
|
! P1> Precipitation evaporation/sublimation |
398 |
|
|
! -------------------------------------------------------------------- |
399 |
|
|
! A part of the precipitation coming from above is evaporated/sublimated. |
400 |
|
|
! -------------------------------------------------------------------- |
401 |
|
|
|
402 |
✓✗ |
11232 |
IF (iflag_evap_prec.GE.1) THEN |
403 |
|
|
|
404 |
|
|
! Calculation of saturation specific humidity |
405 |
|
|
! depending on temperature: |
406 |
|
11232 |
CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:)) |
407 |
|
|
! wrt liquid water |
408 |
|
11232 |
CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,1,.false.,qsl(:),dqsl(:)) |
409 |
|
|
! wrt ice |
410 |
|
11232 |
CALL calc_qsat_ecmwf(klon,zt(:),qzero(:),pplay(:,k),RTT,2,.false.,qsi(:),dqsi(:)) |
411 |
|
|
|
412 |
✓✓ |
11175840 |
DO i = 1, klon |
413 |
|
|
|
414 |
|
|
! if precipitation |
415 |
✓✓ |
11175840 |
IF (zrfl(i)+zifl(i).GT.0.) THEN |
416 |
|
|
|
417 |
|
|
! LudoTP: we only account for precipitation evaporation in the clear-sky (iflag_evap_prec>=4). |
418 |
|
|
! c_iso: likely important to distinguish cs from neb isotope precipitation |
419 |
|
|
|
420 |
✓✗ |
3250315 |
IF (iflag_evap_prec.GE.4) THEN |
421 |
|
3250315 |
zrfl(i) = zrflclr(i) |
422 |
|
3250315 |
zifl(i) = ziflclr(i) |
423 |
|
|
ENDIF |
424 |
|
|
|
425 |
✗✓ |
3250315 |
IF (iflag_evap_prec.EQ.1) THEN |
426 |
|
|
znebprecip(i)=zneb(i) |
427 |
|
|
ELSE |
428 |
|
3250315 |
znebprecip(i)=MAX(zneb(i),znebprecip(i)) |
429 |
|
|
ENDIF |
430 |
|
|
|
431 |
✓✗ |
3250315 |
IF (iflag_evap_prec.GT.4) THEN |
432 |
|
|
! Max evaporation not to saturate the clear sky precip fraction |
433 |
|
|
! i.e. the fraction where evaporation occurs |
434 |
|
3250315 |
zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecipclr(i)) |
435 |
|
|
ELSEIF (iflag_evap_prec .EQ. 4) THEN |
436 |
|
|
! Max evaporation not to saturate the whole mesh |
437 |
|
|
! Pay attention -> lead to unrealistic and excessive evaporation |
438 |
|
|
zqev0 = MAX(0.0, zqs(i)-zq(i)) |
439 |
|
|
ELSE |
440 |
|
|
! Max evap not to saturate the fraction below the cloud |
441 |
|
|
zqev0 = MAX(0.0, (zqs(i)-zq(i))*znebprecip(i)) |
442 |
|
|
ENDIF |
443 |
|
|
|
444 |
|
|
! Evaporation of liquid precipitation coming from above |
445 |
|
|
! dP/dz=beta*(1-q/qsat)*sqrt(P) |
446 |
|
|
! formula from Sundquist 1988, Klemp & Wilhemson 1978 |
447 |
|
|
! LTP: evaporation only in the clear sky part (iflag_evap_prec>=4) |
448 |
|
|
|
449 |
✗✓ |
3250315 |
IF (iflag_evap_prec.EQ.3) THEN |
450 |
|
|
zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl(i)) & |
451 |
|
|
*SQRT(zrfl(i)/max(1.e-4,znebprecip(i))) & |
452 |
|
|
*(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG |
453 |
✓✗ |
3250315 |
ELSE IF (iflag_evap_prec.GE.4) THEN |
454 |
|
|
zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl(i)) & |
455 |
|
|
*SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) & |
456 |
|
3250315 |
*(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG |
457 |
|
|
ELSE |
458 |
|
|
zqevt = 1.*coef_eva*(1.0-zq(i)/qsl(i))*SQRT(zrfl(i)) & |
459 |
|
|
*(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG |
460 |
|
|
ENDIF |
461 |
|
|
|
462 |
|
|
zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) & |
463 |
|
3250315 |
*RG*dtime/(paprs(i,k)-paprs(i,k+1)) |
464 |
|
|
|
465 |
|
|
! sublimation of the solid precipitation coming from above |
466 |
✗✓ |
3250315 |
IF (iflag_evap_prec.EQ.3) THEN |
467 |
|
|
zqevti = znebprecip(i)*coef_eva_i*(1.0-zq(i)/qsi(i)) & |
468 |
|
|
*SQRT(zifl(i)/max(1.e-4,znebprecip(i))) & |
469 |
|
|
*(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG |
470 |
✓✗ |
3250315 |
ELSE IF (iflag_evap_prec.GE.4) THEN |
471 |
|
|
zqevti = znebprecipclr(i)*coef_eva_i*(1.0-zq(i)/qsi(i)) & |
472 |
|
|
*SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) & |
473 |
|
3250315 |
*(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG |
474 |
|
|
ELSE |
475 |
|
|
zqevti = 1.*coef_eva_i*(1.0-zq(i)/qsi(i))*SQRT(zifl(i)) & |
476 |
|
|
*(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG |
477 |
|
|
ENDIF |
478 |
|
|
|
479 |
|
|
zqevti = MAX(0.0,MIN(zqevti,zifl(i))) & |
480 |
|
3250315 |
*RG*dtime/(paprs(i,k)-paprs(i,k+1)) |
481 |
|
|
|
482 |
|
|
! A. JAM |
483 |
|
|
! Evaporation limit: we ensure that the layer's fraction below |
484 |
|
|
! the cloud or the whole mesh (depending on iflag_evap_prec) |
485 |
|
|
! does not reach saturation. In this case, we |
486 |
|
|
! redistribute zqev0 conserving the ratio liquid/ice |
487 |
|
|
|
488 |
|
|
! todoan: check the consistency of this formula |
489 |
✓✓ |
3250315 |
IF (zqevt+zqevti.GT.zqev0) THEN |
490 |
|
156909 |
zqev=zqev0*zqevt/(zqevt+zqevti) |
491 |
|
156909 |
zqevi=zqev0*zqevti/(zqevt+zqevti) |
492 |
|
|
ELSE |
493 |
|
|
zqev=zqevt |
494 |
|
|
zqevi=zqevti |
495 |
|
|
ENDIF |
496 |
|
|
|
497 |
|
|
|
498 |
|
|
! New solid and liquid precipitation fluxes after evap and sublimation |
499 |
|
|
zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) & |
500 |
|
3250315 |
/RG/dtime) |
501 |
|
|
zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) & |
502 |
|
3250315 |
/RG/dtime) |
503 |
|
|
|
504 |
|
|
|
505 |
|
|
! vapor, temperature, precip fluxes update |
506 |
|
|
zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) & |
507 |
|
3250315 |
* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime |
508 |
|
|
zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) & |
509 |
|
3250315 |
* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime |
510 |
|
|
zt(i) = zt(i) + (zrfln(i)-zrfl(i)) & |
511 |
|
|
* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime & |
512 |
|
|
* RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) & |
513 |
|
|
+ (zifln(i)-zifl(i)) & |
514 |
|
|
* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime & |
515 |
|
3250315 |
* RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) |
516 |
|
|
|
517 |
|
|
! New values of liquid and solid precipitation |
518 |
|
3250315 |
zrfl(i) = zrfln(i) |
519 |
|
3250315 |
zifl(i) = zifln(i) |
520 |
|
|
|
521 |
|
|
! c_iso here call_reevap that updates isotopic zrfl, zifl (in inout) |
522 |
|
|
! due to evap + sublim |
523 |
|
|
|
524 |
|
|
|
525 |
✓✗ |
3250315 |
IF (iflag_evap_prec.GE.4) THEN |
526 |
|
3250315 |
zrflclr(i) = zrfl(i) |
527 |
|
3250315 |
ziflclr(i) = zifl(i) |
528 |
✓✓ |
3250315 |
IF(zrflclr(i) + ziflclr(i).LE.0) THEN |
529 |
|
2516867 |
znebprecipclr(i) = 0.0 |
530 |
|
|
ENDIF |
531 |
|
3250315 |
zrfl(i) = zrflclr(i) + zrflcld(i) |
532 |
|
3250315 |
zifl(i) = ziflclr(i) + ziflcld(i) |
533 |
|
|
ENDIF |
534 |
|
|
|
535 |
|
|
! c_iso duplicate for isotopes or loop on isotopes |
536 |
|
|
|
537 |
|
|
! Melting: |
538 |
|
3250315 |
zmelt = ((zt(i)-RTT)/(ztfondue-RTT)) ! JYG |
539 |
|
|
! precip fraction that is melted |
540 |
|
3250315 |
zmelt = MIN(MAX(zmelt,0.),1.) |
541 |
|
|
|
542 |
|
|
! update of rainfall and snowfall due to melting |
543 |
✓✗ |
3250315 |
IF (iflag_evap_prec.GE.4) THEN |
544 |
|
3250315 |
zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i) |
545 |
|
3250315 |
zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i) |
546 |
|
3250315 |
zrfl(i)=zrflclr(i)+zrflcld(i) |
547 |
|
|
|
548 |
|
3250315 |
ziflclr(i)=ziflclr(i)*(1.-zmelt) |
549 |
|
3250315 |
ziflcld(i)=ziflcld(i)*(1.-zmelt) |
550 |
|
3250315 |
zifl(i)=ziflclr(i)+ziflcld(i) |
551 |
|
|
|
552 |
|
|
ELSE |
553 |
|
|
zrfl(i)=zrfl(i)+zmelt*zifl(i) |
554 |
|
|
|
555 |
|
|
zifl(i)=zifl(i)*(1.-zmelt) |
556 |
|
|
ENDIF |
557 |
|
|
|
558 |
|
|
|
559 |
|
|
! c_iso: melting of isotopic precipi with zmelt (no fractionation) |
560 |
|
|
|
561 |
|
|
! Latent heat of melting with precipitation thermalization |
562 |
|
|
zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) & |
563 |
|
3250315 |
*RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) |
564 |
|
|
|
565 |
|
|
|
566 |
|
|
ELSE |
567 |
|
|
! if no precip, we reinitialize the cloud fraction used for the precip to 0 |
568 |
|
7914293 |
znebprecip(i)=0. |
569 |
|
|
|
570 |
|
|
ENDIF ! (zrfl(i)+zifl(i).GT.0.) |
571 |
|
|
|
572 |
|
|
ENDDO |
573 |
|
|
|
574 |
|
|
ENDIF ! (iflag_evap_prec>=1) |
575 |
|
|
|
576 |
|
|
! -------------------------------------------------------------------- |
577 |
|
|
! End precip evaporation |
578 |
|
|
! -------------------------------------------------------------------- |
579 |
|
|
|
580 |
|
|
! Calculation of qsat, L/Cp*dqsat/dT and ncoreczq counter |
581 |
|
|
!------------------------------------------------------- |
582 |
|
|
|
583 |
✓✓ |
11175840 |
qtot(:)=zq(:)+zmqc(:) |
584 |
|
11232 |
CALL calc_qsat_ecmwf(klon,zt(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:)) |
585 |
✓✓ |
11175840 |
DO i = 1, klon |
586 |
|
11164608 |
zdelta = MAX(0.,SIGN(1.,RTT-zt(i))) |
587 |
|
11164608 |
zdqsdT_raw(i) = zdqs(i)*RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta) |
588 |
✗✓ |
11175840 |
IF (zq(i) .LT. 1.e-15) THEN |
589 |
|
|
ncoreczq=ncoreczq+1 |
590 |
|
|
zq(i)=1.e-15 |
591 |
|
|
ENDIF |
592 |
|
|
! c_iso: do something similar for isotopes |
593 |
|
|
|
594 |
|
|
ENDDO |
595 |
|
|
|
596 |
|
|
! -------------------------------------------------------------------- |
597 |
|
|
! P2> Cloud formation |
598 |
|
|
!--------------------------------------------------------------------- |
599 |
|
|
! |
600 |
|
|
! Unlike fisrtilp, we always assume a 'fractional cloud' approach |
601 |
|
|
! i.e. clouds occupy only a fraction of the mesh (the subgrid distribution |
602 |
|
|
! is prescribed and depends on large scale variables and boundary layer |
603 |
|
|
! properties) |
604 |
|
|
! The decrease in condensed part due tu latent heating is taken into |
605 |
|
|
! account |
606 |
|
|
! ------------------------------------------------------------------- |
607 |
|
|
|
608 |
|
|
! P2.1> With the PDFs (log-normal, bigaussian) |
609 |
|
|
! cloud properties calculation with the initial values of t and q |
610 |
|
|
! ---------------------------------------------------------------- |
611 |
|
|
|
612 |
|
|
! initialise gammasat and qincloud_mpc |
613 |
✓✓ |
11175840 |
gammasat(:)=1. |
614 |
✓✓ |
11175840 |
qincloud_mpc(:)=0. |
615 |
|
|
|
616 |
✓✗ |
11232 |
IF (iflag_cld_th.GE.5) THEN |
617 |
|
|
|
618 |
✓✗ |
11232 |
IF (iflag_mpc_bl .LT. 1) THEN |
619 |
|
|
|
620 |
✗✓ |
11232 |
IF (iflag_cloudth_vert.LE.2) THEN |
621 |
|
|
|
622 |
|
|
CALL cloudth(klon,klev,k,ztv, & |
623 |
|
|
zq,zqta,fraca, & |
624 |
|
|
qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, & |
625 |
|
|
ratqs,zqs,t) |
626 |
|
|
|
627 |
✓✗ |
11232 |
ELSEIF (iflag_cloudth_vert.GE.3 .AND. iflag_cloudth_vert.LE.5) THEN |
628 |
|
|
|
629 |
|
|
CALL cloudth_v3(klon,klev,k,ztv, & |
630 |
|
|
zq,zqta,fraca, & |
631 |
|
|
qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & |
632 |
|
11232 |
ratqs,zqs,t) |
633 |
|
|
|
634 |
|
|
!Jean Jouhaud's version, Decembre 2018 |
635 |
|
|
!------------------------------------- |
636 |
|
|
|
637 |
|
|
ELSEIF (iflag_cloudth_vert.EQ.6) THEN |
638 |
|
|
|
639 |
|
|
CALL cloudth_v6(klon,klev,k,ztv, & |
640 |
|
|
zq,zqta,fraca, & |
641 |
|
|
qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & |
642 |
|
|
ratqs,zqs,t) |
643 |
|
|
|
644 |
|
|
ENDIF |
645 |
|
|
|
646 |
|
|
ELSE |
647 |
|
|
! cloudth_v3 + cold microphysical considerations |
648 |
|
|
! + stationary mixed-phase cloud model |
649 |
|
|
|
650 |
|
|
CALL cloudth_mpc(klon,klev,k,iflag_mpc_bl,mpc_bl_points, & |
651 |
|
|
t,ztv,zq,zqta,fraca, zpspsk, & |
652 |
|
|
paprs,pplay,ztla,zthl,ratqs,zqs,psfl, & |
653 |
|
|
qcloud,qincloud_mpc,icefrac_mpc,ctot,ctot_vol) |
654 |
|
|
|
655 |
|
|
ENDIF ! iflag_mpc_bl |
656 |
|
|
|
657 |
✓✓ |
11175840 |
DO i=1,klon |
658 |
|
11164608 |
rneb(i,k)=ctot(i,k) |
659 |
|
11164608 |
rneblsvol(i,k)=ctot_vol(i,k) |
660 |
|
11175840 |
zqn(i)=qcloud(i) |
661 |
|
|
ENDDO |
662 |
|
|
|
663 |
|
|
ENDIF |
664 |
|
|
|
665 |
✗✓ |
11232 |
IF (iflag_cld_th .LE. 4) THEN |
666 |
|
|
|
667 |
|
|
! lognormal |
668 |
|
|
lognormale = .TRUE. |
669 |
|
|
|
670 |
✓✗ |
11232 |
ELSEIF (iflag_cld_th .GE. 6) THEN |
671 |
|
|
|
672 |
|
|
! lognormal distribution when no thermals |
673 |
✓✓ |
11175840 |
lognormale = fraca(:,k) < 1e-10 |
674 |
|
|
|
675 |
|
|
ELSE |
676 |
|
|
! When iflag_cld_th=5, we always assume |
677 |
|
|
! bi-gaussian distribution |
678 |
|
|
lognormale = .FALSE. |
679 |
|
|
|
680 |
|
|
ENDIF |
681 |
|
|
|
682 |
✓✓ |
11175840 |
DT(:) = 0. |
683 |
✓✓ |
11175840 |
n_i(:)=0 |
684 |
✓✓ |
11175840 |
Tbef(:)=zt(:) |
685 |
✓✓ |
11175840 |
qlbef(:)=0. |
686 |
|
|
|
687 |
|
|
! Treatment of non-boundary layer clouds (lognormale) |
688 |
|
|
! condensation with qsat(T) variation (adaptation) |
689 |
|
|
! Iterative resolution to converge towards qsat |
690 |
|
|
! with update of temperature, ice fraction and qsat at |
691 |
|
|
! each iteration |
692 |
|
|
|
693 |
|
|
! todoan -> sensitivity to iflag_fisrtilp_qsat |
694 |
✓✓ |
67392 |
DO iter=1,iflag_fisrtilp_qsat+1 |
695 |
|
|
|
696 |
✓✓ |
55879200 |
DO i=1,klon |
697 |
|
|
|
698 |
|
|
! keepgoing = .true. while convergence is not satisfied |
699 |
|
55823040 |
keepgoing(i)=ABS(DT(i)).GT.DDT0 |
700 |
|
|
|
701 |
✓✓✓✓
✓✓ |
55879200 |
IF ((keepgoing(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN |
702 |
|
|
|
703 |
|
|
! if not convergence: |
704 |
|
|
|
705 |
|
|
! P2.2.1> cloud fraction and condensed water mass calculation |
706 |
|
|
! Calculated variables: |
707 |
|
|
! rneb : cloud fraction |
708 |
|
|
! zqn : total water within the cloud |
709 |
|
|
! zcond: mean condensed water within the mesh |
710 |
|
|
! rhcl: clear-sky relative humidity |
711 |
|
|
!--------------------------------------------------------------- |
712 |
|
|
|
713 |
|
|
! new temperature that only serves in the iteration process: |
714 |
|
11599051 |
Tbef(i)=Tbef(i)+DT(i) |
715 |
|
|
|
716 |
|
|
! Rneb, qzn and zcond for lognormal PDFs |
717 |
|
11599051 |
qtot(i)=zq(i)+zmqc(i) |
718 |
|
|
|
719 |
|
|
ENDIF |
720 |
|
|
|
721 |
|
|
ENDDO |
722 |
|
|
|
723 |
|
|
! Calculation of saturation specific humidity and ice fraction |
724 |
|
56160 |
CALL calc_qsat_ecmwf(klon,Tbef(:),qtot(:),pplay(:,k),RTT,0,.false.,zqs(:),zdqs(:)) |
725 |
|
56160 |
CALL calc_gammasat(klon,Tbef(:),qtot(:),pplay(:,k),gammasat(:),dgammasatdt(:)) |
726 |
|
|
! saturation may occur at a humidity different from qsat (gamma qsat), so gamma correction for dqs |
727 |
✓✓ |
55879200 |
zdqs(:) = gammasat(:)*zdqs(:)+zqs(:)*dgammasatdt(:) |
728 |
|
|
! cloud phase determination |
729 |
✗✓ |
56160 |
IF (iflag_t_glace.GE.4) THEN |
730 |
|
|
! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top |
731 |
|
|
CALL distance_to_cloud_top(klon,klev,k,t,pplay,paprs,rneb,distcltop1D) |
732 |
|
|
ENDIF |
733 |
|
56160 |
CALL icefrac_lscp(klon, zt(:), iflag_ice_thermo, distcltop1D(:),zfice(:),dzfice(:)) |
734 |
|
|
|
735 |
✓✓ |
55890432 |
DO i=1,klon !todoan : check if loop in i is needed |
736 |
|
|
|
737 |
✓✓✓✓
✓✓ |
55879200 |
IF ((keepgoing(i) .OR. (n_i(i) .EQ. 0)) .AND. lognormale(i)) THEN |
738 |
|
|
|
739 |
|
11599051 |
zpdf_sig(i)=ratqs(i,k)*zq(i) |
740 |
|
11599051 |
zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2)) |
741 |
|
11599051 |
zpdf_delta(i)=log(zq(i)/(gammasat(i)*zqs(i))) |
742 |
|
11599051 |
zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.)) |
743 |
|
11599051 |
zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.)) |
744 |
|
11599051 |
zpdf_e1(i)=zpdf_a(i)-zpdf_b(i) |
745 |
|
11599051 |
zpdf_e1(i)=sign(min(ABS(zpdf_e1(i)),5.),zpdf_e1(i)) |
746 |
|
11599051 |
zpdf_e1(i)=1.-erf(zpdf_e1(i)) |
747 |
|
11599051 |
zpdf_e2(i)=zpdf_a(i)+zpdf_b(i) |
748 |
|
11599051 |
zpdf_e2(i)=sign(min(ABS(zpdf_e2(i)),5.),zpdf_e2(i)) |
749 |
|
11599051 |
zpdf_e2(i)=1.-erf(zpdf_e2(i)) |
750 |
|
|
|
751 |
✗✓✗✗
|
11599051 |
IF ((.NOT.ok_ice_sursat).OR.(Tbef(i).GT.t_glace_min)) THEN |
752 |
|
|
|
753 |
✓✓ |
11599051 |
IF (zpdf_e1(i).LT.1.e-10) THEN |
754 |
|
7219503 |
rneb(i,k)=0. |
755 |
|
7219503 |
zqn(i)=gammasat(i)*zqs(i) |
756 |
|
|
ELSE |
757 |
|
4379548 |
rneb(i,k)=0.5*zpdf_e1(i) |
758 |
|
4379548 |
zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i) |
759 |
|
|
ENDIF |
760 |
|
|
|
761 |
|
11599051 |
rnebss(i,k)=0.0 !--added by OB (needed because of the convergence loop on time) |
762 |
|
11599051 |
fcontrN(i,k)=0.0 !--idem |
763 |
|
11599051 |
fcontrP(i,k)=0.0 !--idem |
764 |
|
11599051 |
qss(i,k)=0.0 !--idem |
765 |
|
|
|
766 |
|
|
ELSE ! in case of ice supersaturation by Audran |
767 |
|
|
|
768 |
|
|
!------------------------------------ |
769 |
|
|
! ICE SUPERSATURATION |
770 |
|
|
!------------------------------------ |
771 |
|
|
|
772 |
|
|
CALL ice_sursat(pplay(i,k), paprs(i,k)-paprs(i,k+1), dtime, i, k, t(i,k), zq(i), & |
773 |
|
|
gamma_ss(i,k), zqs(i), Tbef(i), rneb_seri(i,k), ratqs(i,k), & |
774 |
|
|
rneb(i,k), zqn(i), rnebss(i,k), qss(i,k), & |
775 |
|
|
Tcontr(i,k), qcontr(i,k), qcontr2(i,k), fcontrN(i,k), fcontrP(i,k) ) |
776 |
|
|
|
777 |
|
|
ENDIF ! ((flag_ice_sursat.eq.0).or.(Tbef(i).gt.t_glace_min)) |
778 |
|
|
|
779 |
|
|
! If vertical heterogeneity, change fraction by volume as well |
780 |
✓✗ |
11599051 |
IF (iflag_cloudth_vert.GE.3) THEN |
781 |
|
11599051 |
ctot_vol(i,k)=rneb(i,k) |
782 |
|
11599051 |
rneblsvol(i,k)=ctot_vol(i,k) |
783 |
|
|
ENDIF |
784 |
|
|
|
785 |
|
|
|
786 |
|
|
! P2.2.2> Approximative calculation of temperature variation DT |
787 |
|
|
! due to condensation. |
788 |
|
|
! Calculated variables: |
789 |
|
|
! dT : temperature change due to condensation |
790 |
|
|
!--------------------------------------------------------------- |
791 |
|
|
|
792 |
|
|
|
793 |
✓✓ |
11599051 |
IF (zfice(i).LT.1) THEN |
794 |
|
4401650 |
cste=RLVTT |
795 |
|
|
ELSE |
796 |
|
7197401 |
cste=RLSTT |
797 |
|
|
ENDIF |
798 |
|
|
|
799 |
|
11599051 |
qlbef(i)=max(0.,zqn(i)-zqs(i)) |
800 |
|
|
num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT & |
801 |
|
11599051 |
+zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i) |
802 |
|
|
denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) & |
803 |
|
|
-(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k) & |
804 |
|
11599051 |
*qlbef(i)*dzfice(i) |
805 |
|
|
! here we update a provisory temperature variable that only serves in the iteration |
806 |
|
|
! process |
807 |
|
11599051 |
DT(i)=num/denom |
808 |
|
11599051 |
n_i(i)=n_i(i)+1 |
809 |
|
|
|
810 |
|
|
ENDIF ! end keepgoing |
811 |
|
|
|
812 |
|
|
ENDDO ! end loop on i |
813 |
|
|
|
814 |
|
|
ENDDO ! iter=1,iflag_fisrtilp_qsat+1 |
815 |
|
|
|
816 |
|
|
! P2.3> Final quantities calculation |
817 |
|
|
! Calculated variables: |
818 |
|
|
! rneb : cloud fraction |
819 |
|
|
! zcond: mean condensed water in the mesh |
820 |
|
|
! zqn : mean water vapor in the mesh |
821 |
|
|
! zfice: ice fraction in clouds |
822 |
|
|
! zt : temperature |
823 |
|
|
! rhcl : clear-sky relative humidity |
824 |
|
|
! ---------------------------------------------------------------- |
825 |
|
|
|
826 |
|
|
|
827 |
|
|
! For iflag_t_glace GE 4 the phase partition function dependends on temperature AND distance to cloud top |
828 |
✗✓ |
11232 |
IF (iflag_t_glace.GE.4) THEN |
829 |
|
|
CALL distance_to_cloud_top(klon,klev,k,t,pplay,paprs,rneb,distcltop1D) |
830 |
|
|
distcltop(:,k)=distcltop1D(:) |
831 |
|
|
ENDIF |
832 |
|
|
|
833 |
|
|
! Partition function in stratiform clouds (will be overwritten in boundary-layer MPCs) |
834 |
|
11232 |
CALL icefrac_lscp(klon,zt(:),iflag_ice_thermo,distcltop1D(:),zfice(:), dzfice(:)) |
835 |
|
|
|
836 |
|
|
|
837 |
|
|
! Water vapor update, Phase determination and subsequent latent heat exchange |
838 |
✓✓ |
11175840 |
DO i=1, klon |
839 |
|
|
|
840 |
✗✓ |
11164608 |
IF (mpc_bl_points(i,k) .GT. 0) THEN |
841 |
|
|
zcond(i) = MAX(0.0,qincloud_mpc(i))*rneb(i,k) |
842 |
|
|
! following line is very strange and probably wrong |
843 |
|
|
rhcl(i,k)= (zqs(i)+zq(i))/2./zqs(i) |
844 |
|
|
! water vapor update and partition function if thermals |
845 |
|
|
zq(i) = zq(i) - zcond(i) |
846 |
|
|
zfice(i)=icefrac_mpc(i,k) |
847 |
|
|
|
848 |
|
|
ELSE |
849 |
|
|
|
850 |
|
|
! Checks on rneb, rhcl and zqn |
851 |
✓✓ |
11164608 |
IF (rneb(i,k) .LE. 0.0) THEN |
852 |
|
7675015 |
zqn(i) = 0.0 |
853 |
|
7675015 |
rneb(i,k) = 0.0 |
854 |
|
7675015 |
zcond(i) = 0.0 |
855 |
|
7675015 |
rhcl(i,k)=zq(i)/zqs(i) |
856 |
✓✓ |
3489593 |
ELSE IF (rneb(i,k) .GE. 1.0) THEN |
857 |
|
16931 |
zqn(i) = zq(i) |
858 |
|
16931 |
rneb(i,k) = 1.0 |
859 |
|
16931 |
zcond(i) = MAX(0.0,zqn(i)-gammasat(i)*zqs(i)) |
860 |
|
16931 |
rhcl(i,k)=1.0 |
861 |
|
|
ELSE |
862 |
|
3472662 |
zcond(i) = MAX(0.0,zqn(i)-gammasat(i)*zqs(i))*rneb(i,k) |
863 |
|
|
! following line is very strange and probably wrong: |
864 |
|
3472662 |
rhcl(i,k)=(zqs(i)+zq(i))/2./zqs(i) |
865 |
|
|
ENDIF |
866 |
|
|
|
867 |
|
|
! water vapor update |
868 |
|
11164608 |
zq(i) = zq(i) - zcond(i) |
869 |
|
|
|
870 |
|
|
ENDIF |
871 |
|
|
|
872 |
|
|
! c_iso : routine that computes in-cloud supersaturation |
873 |
|
|
! c_iso condensation of isotopes (zcond, zsursat, zfice, zq in input) |
874 |
|
|
|
875 |
|
|
! temperature update due to phase change |
876 |
|
|
zt(i) = zt(i) + (1.-zfice(i))*zcond(i) & |
877 |
|
|
& * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) & |
878 |
|
11175840 |
+zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) |
879 |
|
|
ENDDO |
880 |
|
|
|
881 |
|
|
! If vertical heterogeneity, change volume fraction |
882 |
✓✗ |
11232 |
IF (iflag_cloudth_vert .GE. 3) THEN |
883 |
✓✓ |
11175840 |
ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.) |
884 |
✓✓ |
11175840 |
rneblsvol(1:klon,k)=ctot_vol(1:klon,k) |
885 |
|
|
ENDIF |
886 |
|
|
|
887 |
|
|
! ---------------------------------------------------------------- |
888 |
|
|
! P3> Precipitation formation |
889 |
|
|
! ---------------------------------------------------------------- |
890 |
|
|
|
891 |
|
|
! LTP: |
892 |
✓✗ |
11232 |
IF (iflag_evap_prec .GE. 4) THEN |
893 |
|
|
|
894 |
|
|
!Partitionning between precipitation coming from clouds and that coming from CS |
895 |
|
|
|
896 |
|
|
!0) Calculate tot_zneb, total cloud fraction above the cloud |
897 |
|
|
!assuming a maximum-random overlap (voir Jakob and Klein, 2000) |
898 |
|
|
|
899 |
✓✓ |
11175840 |
DO i=1, klon |
900 |
|
|
tot_znebn(i) = 1. - (1.-tot_zneb(i))*(1 - max(rneb(i,k),zneb(i))) & |
901 |
|
11164608 |
/(1.-min(zneb(i),1.-smallestreal)) |
902 |
|
11164608 |
d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i) |
903 |
|
11164608 |
tot_zneb(i) = tot_znebn(i) |
904 |
|
|
|
905 |
|
|
|
906 |
|
|
!1) Cloudy to clear air |
907 |
|
11164608 |
d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i,k),znebprecipcld(i)) |
908 |
✓✓ |
11164608 |
IF (znebprecipcld(i) .GT. 0.) THEN |
909 |
|
2677825 |
d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i) |
910 |
|
2677825 |
d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i) |
911 |
|
|
ELSE |
912 |
|
8486783 |
d_zrfl_cld_clr(i) = 0. |
913 |
|
8486783 |
d_zifl_cld_clr(i) = 0. |
914 |
|
|
ENDIF |
915 |
|
|
|
916 |
|
|
!2) Clear to cloudy air |
917 |
|
11164608 |
d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i,k) - d_tot_zneb(i) - zneb(i))) |
918 |
✓✓ |
11164608 |
IF (znebprecipclr(i) .GT. 0) THEN |
919 |
|
733448 |
d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i) |
920 |
|
733448 |
d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i) |
921 |
|
|
ELSE |
922 |
|
10431160 |
d_zrfl_clr_cld(i) = 0. |
923 |
|
10431160 |
d_zifl_clr_cld(i) = 0. |
924 |
|
|
ENDIF |
925 |
|
|
|
926 |
|
|
!Update variables |
927 |
|
11164608 |
znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i) |
928 |
|
11164608 |
znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i) |
929 |
|
11164608 |
zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i) |
930 |
|
11164608 |
ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i) |
931 |
|
11164608 |
zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i) |
932 |
|
11175840 |
ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i) |
933 |
|
|
|
934 |
|
|
! c_iso do the same thing for isotopes precip |
935 |
|
|
ENDDO |
936 |
|
|
ENDIF |
937 |
|
|
|
938 |
|
|
|
939 |
|
|
! Autoconversion |
940 |
|
|
! ------------------------------------------------------------------------------- |
941 |
|
|
|
942 |
|
|
|
943 |
|
|
! Initialisation of zoliq and radocond variables |
944 |
|
|
|
945 |
✓✓ |
11175840 |
DO i = 1, klon |
946 |
|
11164608 |
zoliq(i) = zcond(i) |
947 |
|
11164608 |
zoliqi(i)= zoliq(i)*zfice(i) |
948 |
|
11164608 |
zoliql(i)= zoliq(i)*(1.-zfice(i)) |
949 |
|
|
! c_iso : initialisation of zoliq* also for isotopes |
950 |
|
11164608 |
zrho(i,k) = pplay(i,k) / zt(i) / RD |
951 |
|
11164608 |
zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i,k)*RG) |
952 |
|
11164608 |
iwc(i) = 0. |
953 |
|
11164608 |
zneb(i) = MAX(rneb(i,k),seuil_neb) |
954 |
|
11164608 |
radocond(i,k) = zoliq(i)/REAL(niter_lscp+1) |
955 |
|
11164608 |
radicefrac(i,k) = zfice(i) |
956 |
|
11164608 |
radocondi(i,k)=zoliq(i)*zfice(i)/REAL(niter_lscp+1) |
957 |
|
11175840 |
radocondl(i,k)=zoliq(i)*(1.-zfice(i))/REAL(niter_lscp+1) |
958 |
|
|
ENDDO |
959 |
|
|
|
960 |
|
|
|
961 |
✓✓ |
67392 |
DO n = 1, niter_lscp |
962 |
|
|
|
963 |
✓✓ |
55879200 |
DO i=1,klon |
964 |
✓✓ |
55879200 |
IF (rneb(i,k).GT.0.0) THEN |
965 |
|
17447965 |
iwc(i) = zrho(i,k) * zoliqi(i) / zneb(i) ! in-cloud ice condensate content |
966 |
|
|
ENDIF |
967 |
|
|
ENDDO |
968 |
|
|
|
969 |
|
56160 |
CALL fallice_velocity(klon,iwc(:),zt(:),zrho(:,k),pplay(:,k),ptconv(:,k),velo(:,k)) |
970 |
|
|
|
971 |
✓✓ |
55890432 |
DO i = 1, klon |
972 |
|
|
|
973 |
✓✓ |
55879200 |
IF (rneb(i,k).GT.0.0) THEN |
974 |
|
|
|
975 |
|
|
! Initialization of zrain, zsnow and zprecip: |
976 |
|
|
zrain=0. |
977 |
|
|
zsnow=0. |
978 |
|
|
zprecip=0. |
979 |
|
|
! c_iso same init for isotopes. Externalisation? |
980 |
|
|
|
981 |
✓✓ |
17447965 |
IF (zneb(i).EQ.seuil_neb) THEN |
982 |
|
|
zprecip = 0.0 |
983 |
|
|
zsnow = 0.0 |
984 |
|
|
zrain= 0.0 |
985 |
|
|
ELSE |
986 |
|
|
|
987 |
✓✓ |
12022170 |
IF (ptconv(i,k)) THEN ! if convective point |
988 |
|
1223420 |
zcl=cld_lc_con |
989 |
|
1223420 |
zct=1./cld_tau_con |
990 |
|
1223420 |
zexpo=cld_expo_con |
991 |
|
1223420 |
ffallv=ffallv_con |
992 |
|
|
ELSE |
993 |
|
10798750 |
zcl=cld_lc_lsc |
994 |
|
10798750 |
zct=1./cld_tau_lsc |
995 |
|
10798750 |
zexpo=cld_expo_lsc |
996 |
|
10798750 |
ffallv=ffallv_lsc |
997 |
|
|
ENDIF |
998 |
|
|
|
999 |
|
|
|
1000 |
|
|
! if vertical heterogeneity is taken into account, we use |
1001 |
|
|
! the "true" volume fraction instead of a modified |
1002 |
|
|
! surface fraction (which is larger and artificially |
1003 |
|
|
! reduces the in-cloud water). |
1004 |
|
|
|
1005 |
|
|
! Liquid water quantity to remove: zchau (Sundqvist, 1978) |
1006 |
|
|
! dqliq/dt=-qliq/tau*(1-exp(-qcin/clw)**2) |
1007 |
|
|
!......................................................... |
1008 |
✓✗✗✓
|
12022170 |
IF ((iflag_cloudth_vert.GE.3).AND.(iflag_rain_incloud_vol.EQ.1)) THEN |
1009 |
|
|
|
1010 |
|
|
! if vertical heterogeneity is taken into account, we use |
1011 |
|
|
! the "true" volume fraction instead of a modified |
1012 |
|
|
! surface fraction (which is larger and artificially |
1013 |
|
|
! reduces the in-cloud water). |
1014 |
|
|
effective_zneb=ctot_vol(i,k) |
1015 |
|
|
ELSE |
1016 |
|
|
effective_zneb=zneb(i) |
1017 |
|
|
ENDIF |
1018 |
|
|
|
1019 |
|
|
|
1020 |
✗✓ |
12022170 |
IF (iflag_autoconversion .EQ. 2) THEN |
1021 |
|
|
! two-steps resolution with niter_lscp=1 sufficient |
1022 |
|
|
! we first treat the second term (with exponential) in an explicit way |
1023 |
|
|
! and then treat the first term (-q/tau) in an exact way |
1024 |
|
|
zchau=zoliql(i)*(1.-exp(-dtime/REAL(niter_lscp)*zct & |
1025 |
|
|
*(1.-exp(-(zoliql(i)/effective_zneb/zcl)**zexpo)))) |
1026 |
|
|
ELSE |
1027 |
|
|
! old explicit resolution with subtimesteps |
1028 |
|
|
zchau = zct*dtime/REAL(niter_lscp)*zoliql(i) & |
1029 |
|
12022170 |
*(1.0-EXP(-(zoliql(i)/effective_zneb/zcl)**zexpo)) |
1030 |
|
|
ENDIF |
1031 |
|
|
|
1032 |
|
|
|
1033 |
|
|
! Ice water quantity to remove (Zender & Kiehl, 1997) |
1034 |
|
|
! dqice/dt=1/rho*d(rho*wice*qice)/dz |
1035 |
|
|
!.................................... |
1036 |
✗✓ |
12022170 |
IF (iflag_autoconversion .EQ. 2) THEN |
1037 |
|
|
! exact resolution, niter_lscp=1 is sufficient but works only |
1038 |
|
|
! with iflag_vice=0 |
1039 |
|
|
IF (zoliqi(i) .GT. 0.) THEN |
1040 |
|
|
zfroi=(zoliqi(i)-((zoliqi(i)**(-dice_velo)) & |
1041 |
|
|
+dice_velo*dtime/REAL(niter_lscp)*cice_velo/zdz(i)*ffallv)**(-1./dice_velo)) |
1042 |
|
|
ELSE |
1043 |
|
|
zfroi=0. |
1044 |
|
|
ENDIF |
1045 |
|
|
ELSE |
1046 |
|
|
! old explicit resolution with subtimesteps |
1047 |
|
12022170 |
zfroi = dtime/REAL(niter_lscp)/zdz(i)*zoliqi(i)*velo(i,k) |
1048 |
|
|
ENDIF |
1049 |
|
|
|
1050 |
|
12022170 |
zrain = MIN(MAX(zchau,0.0),zoliql(i)) |
1051 |
|
12022170 |
zsnow = MIN(MAX(zfroi,0.0),zoliqi(i)) |
1052 |
|
12022170 |
zprecip = MAX(zrain + zsnow,0.0) |
1053 |
|
|
|
1054 |
|
|
ENDIF |
1055 |
|
|
|
1056 |
✓✗ |
17447965 |
IF (iflag_autoconversion .GE. 1) THEN |
1057 |
|
|
! debugged version with phase conservation through the autoconversion process |
1058 |
|
17447965 |
zoliql(i) = MAX(zoliql(i)-1.*zrain , 0.0) |
1059 |
|
17447965 |
zoliqi(i) = MAX(zoliqi(i)-1.*zsnow , 0.0) |
1060 |
|
17447965 |
zoliq(i) = MAX(zoliq(i)-zprecip , 0.0) |
1061 |
|
|
ELSE |
1062 |
|
|
! bugged version with phase resetting |
1063 |
|
|
zoliql(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zrain , 0.0) |
1064 |
|
|
zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zsnow , 0.0) |
1065 |
|
|
zoliq(i) = MAX(zoliq(i)-zprecip , 0.0) |
1066 |
|
|
ENDIF |
1067 |
|
|
|
1068 |
|
|
! c_iso: call isotope_conversion (for readibility) or duplicate |
1069 |
|
|
|
1070 |
|
17447965 |
radocond(i,k) = radocond(i,k) + zoliq(i)/REAL(niter_lscp+1) |
1071 |
|
17447965 |
radocondl(i,k) = radocondl(i,k) + zoliql(i)/REAL(niter_lscp+1) |
1072 |
|
17447965 |
radocondi(i,k) = radocondi(i,k) + zoliqi(i)/REAL(niter_lscp+1) |
1073 |
|
|
|
1074 |
|
|
ENDIF ! rneb >0 |
1075 |
|
|
|
1076 |
|
|
ENDDO ! i = 1,klon |
1077 |
|
|
|
1078 |
|
|
ENDDO ! n = 1,niter |
1079 |
|
|
|
1080 |
|
|
! Precipitation flux calculation |
1081 |
|
|
|
1082 |
✓✓ |
11175840 |
DO i = 1, klon |
1083 |
|
|
|
1084 |
✓✗ |
11164608 |
IF (iflag_evap_prec.GE.4) THEN |
1085 |
|
11164608 |
ziflprev(i)=ziflcld(i) |
1086 |
|
|
ELSE |
1087 |
|
|
ziflprev(i)=zifl(i)*zneb(i) |
1088 |
|
|
ENDIF |
1089 |
|
|
|
1090 |
✓✓ |
11175840 |
IF (rneb(i,k) .GT. 0.0) THEN |
1091 |
|
|
|
1092 |
|
|
! CR&JYG: We account for the Wegener-Findeisen-Bergeron process in the precipitation flux: |
1093 |
|
|
! If T<0C, liquid precip are converted into ice, which leads to an increase in |
1094 |
|
|
! temperature DeltaT. The effect of DeltaT on condensates and precipitation is roughly |
1095 |
|
|
! taken into account through a linearization of the equations and by approximating |
1096 |
|
|
! the liquid precipitation process with a "threshold" process. We assume that |
1097 |
|
|
! condensates are not modified during this operation. Liquid precipitation is |
1098 |
|
|
! removed (in the limit DeltaT<273.15-T). Solid precipitation increases. |
1099 |
|
|
! Water vapor increases as well |
1100 |
|
|
! Note that compared to fisrtilp, we always assume iflag_bergeron=2 |
1101 |
|
|
|
1102 |
|
3489593 |
zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i) |
1103 |
|
3489593 |
zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i)) |
1104 |
|
3489593 |
zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) |
1105 |
|
3489593 |
coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i) |
1106 |
|
|
! Computation of DT if all the liquid precip freezes |
1107 |
|
3489593 |
DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1)) |
1108 |
|
|
! T should not exceed the freezing point |
1109 |
|
|
! that is Delta > RTT-zt(i) |
1110 |
|
3489593 |
DeltaT = max( min( RTT-zt(i), DeltaT) , 0. ) |
1111 |
|
3489593 |
zt(i) = zt(i) + DeltaT |
1112 |
|
|
! water vaporization due to temp. increase |
1113 |
|
3489593 |
Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT |
1114 |
|
|
! we add this vaporized water to the total vapor and we remove it from the precip |
1115 |
|
3489593 |
zq(i) = zq(i) + Deltaq |
1116 |
|
|
! The three "max" lines herebelow protect from rounding errors |
1117 |
|
3489593 |
zcond(i) = max( zcond(i)- Deltaq, 0. ) |
1118 |
|
|
! liquid precipitation converted to ice precip |
1119 |
|
3489593 |
Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT |
1120 |
|
3489593 |
zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. ) |
1121 |
|
|
! iced water budget |
1122 |
|
3489593 |
zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.) |
1123 |
|
|
|
1124 |
|
|
! c_iso : mv here condensation of isotopes + redispatchage en precip |
1125 |
|
|
|
1126 |
✓✗ |
3489593 |
IF (iflag_evap_prec.GE.4) THEN |
1127 |
|
|
zrflcld(i) = zrflcld(i)+zqprecl(i) & |
1128 |
|
3489593 |
*(paprs(i,k)-paprs(i,k+1))/(RG*dtime) |
1129 |
|
|
ziflcld(i) = ziflcld(i)+ zqpreci(i) & |
1130 |
|
3489593 |
*(paprs(i,k)-paprs(i,k+1))/(RG*dtime) |
1131 |
|
3489593 |
znebprecipcld(i) = rneb(i,k) |
1132 |
|
3489593 |
zrfl(i) = zrflcld(i) + zrflclr(i) |
1133 |
|
3489593 |
zifl(i) = ziflcld(i) + ziflclr(i) |
1134 |
|
|
ELSE |
1135 |
|
|
zrfl(i) = zrfl(i)+ zqprecl(i) & |
1136 |
|
|
*(paprs(i,k)-paprs(i,k+1))/(RG*dtime) |
1137 |
|
|
zifl(i) = zifl(i)+ zqpreci(i) & |
1138 |
|
|
*(paprs(i,k)-paprs(i,k+1))/(RG*dtime) |
1139 |
|
|
ENDIF |
1140 |
|
|
! c_iso : same for isotopes |
1141 |
|
|
|
1142 |
|
|
ENDIF ! rneb>0 |
1143 |
|
|
|
1144 |
|
|
ENDDO |
1145 |
|
|
|
1146 |
|
|
! LTP: limit of surface cloud fraction covered by precipitation when the local intensity of the flux is below rain_int_min |
1147 |
|
|
! if iflag_evap_prec>=4 |
1148 |
✓✗ |
11232 |
IF (iflag_evap_prec.GE.4) THEN |
1149 |
|
|
|
1150 |
✓✓ |
11175840 |
DO i=1,klon |
1151 |
|
|
|
1152 |
✓✓ |
11164608 |
IF ((zrflclr(i) + ziflclr(i)) .GT. 0. ) THEN |
1153 |
|
|
znebprecipclr(i) = min(znebprecipclr(i),max(zrflclr(i)/ & |
1154 |
|
1562914 |
(MAX(znebprecipclr(i),seuil_neb)*rain_int_min), ziflclr(i)/(MAX(znebprecipclr(i),seuil_neb)*rain_int_min))) |
1155 |
|
|
ELSE |
1156 |
|
9601694 |
znebprecipclr(i)=0.0 |
1157 |
|
|
ENDIF |
1158 |
|
|
|
1159 |
✓✓ |
11175840 |
IF ((zrflcld(i) + ziflcld(i)) .GT. 0.) THEN |
1160 |
|
|
znebprecipcld(i) = min(znebprecipcld(i), max(zrflcld(i)/ & |
1161 |
|
2777435 |
(MAX(znebprecipcld(i),seuil_neb)*rain_int_min), ziflcld(i)/(MAX(znebprecipcld(i),seuil_neb)*rain_int_min))) |
1162 |
|
|
ELSE |
1163 |
|
8387173 |
znebprecipcld(i)=0.0 |
1164 |
|
|
ENDIF |
1165 |
|
|
|
1166 |
|
|
ENDDO |
1167 |
|
|
|
1168 |
|
|
|
1169 |
|
|
ENDIF |
1170 |
|
|
|
1171 |
|
|
! End of precipitation formation |
1172 |
|
|
! -------------------------------- |
1173 |
|
|
|
1174 |
|
|
! Outputs: |
1175 |
|
|
! Precipitation fluxes at layer interfaces |
1176 |
|
|
! + precipitation fractions + |
1177 |
|
|
! temperature and water species tendencies |
1178 |
✓✓ |
11175840 |
DO i = 1, klon |
1179 |
|
11164608 |
psfl(i,k)=zifl(i) |
1180 |
|
11164608 |
prfl(i,k)=zrfl(i) |
1181 |
|
11164608 |
pfraclr(i,k)=znebprecipclr(i) |
1182 |
|
11164608 |
pfracld(i,k)=znebprecipcld(i) |
1183 |
|
11164608 |
d_ql(i,k) = (1-zfice(i))*zoliq(i) |
1184 |
|
11164608 |
d_qi(i,k) = zfice(i)*zoliq(i) |
1185 |
|
11164608 |
d_q(i,k) = zq(i) - q(i,k) |
1186 |
|
|
! c_iso: same for isotopes |
1187 |
|
11175840 |
d_t(i,k) = zt(i) - t(i,k) |
1188 |
|
|
ENDDO |
1189 |
|
|
|
1190 |
|
|
! Calculation of the concentration of condensates seen by the radiation scheme |
1191 |
|
|
! for liquid phase, we take radocondl |
1192 |
|
|
! for ice phase, we take radocondi if we neglect snowfall, otherwise (ok_radocond_snow=true) |
1193 |
|
|
! we recaulate radocondi to account for contributions from the precipitation flux |
1194 |
|
|
|
1195 |
✗✓✗✗
|
11232 |
IF ((ok_radocond_snow) .AND. (k .LT. klev)) THEN |
1196 |
|
|
! for the solid phase (crystals + snowflakes) |
1197 |
|
|
! we recalculate radocondi by summing |
1198 |
|
|
! the ice content calculated in the mesh |
1199 |
|
|
! + the contribution of the non-evaporated snowfall |
1200 |
|
|
! from the overlying layer |
1201 |
|
|
DO i=1,klon |
1202 |
|
|
IF (ziflprev(i) .NE. 0.0) THEN |
1203 |
|
|
radocondi(i,k)=zoliq(i)*zfice(i)+zqpreci(i)+ziflprev(i)/zrho(i,k+1)/velo(i,k+1) |
1204 |
|
|
ELSE |
1205 |
|
|
radocondi(i,k)=zoliq(i)*zfice(i)+zqpreci(i) |
1206 |
|
|
ENDIF |
1207 |
|
|
radocond(i,k)=radocondl(i,k)+radocondi(i,k) |
1208 |
|
|
ENDDO |
1209 |
|
|
ENDIF |
1210 |
|
|
|
1211 |
|
|
! caculate the percentage of ice in "radocond" so cloud+precip seen by the radiation scheme |
1212 |
✓✓ |
11175840 |
DO i=1,klon |
1213 |
✓✓ |
11175840 |
IF (radocond(i,k) .GT. 0.) THEN |
1214 |
|
3489593 |
radicefrac(i,k)=MIN(MAX(radocondi(i,k)/radocond(i,k),0.),1.) |
1215 |
|
|
ENDIF |
1216 |
|
|
ENDDO |
1217 |
|
|
|
1218 |
|
|
! ---------------------------------------------------------------- |
1219 |
|
|
! P4> Wet scavenging |
1220 |
|
|
! ---------------------------------------------------------------- |
1221 |
|
|
|
1222 |
|
|
!Scavenging through nucleation in the layer |
1223 |
|
|
|
1224 |
✓✓ |
11175840 |
DO i = 1,klon |
1225 |
|
|
|
1226 |
✓✓ |
11164608 |
IF(zcond(i).GT.zoliq(i)+1.e-10) THEN |
1227 |
|
2339358 |
beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime |
1228 |
|
|
ELSE |
1229 |
|
8825250 |
beta(i,k) = 0. |
1230 |
|
|
ENDIF |
1231 |
|
|
|
1232 |
|
11164608 |
zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0)*(paprs(i,k)-paprs(i,k+1))/RG |
1233 |
|
|
|
1234 |
✓✓✓✓
|
11175840 |
IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN |
1235 |
|
|
|
1236 |
✓✓ |
2404434 |
IF (t(i,k) .GE. t_glace_min) THEN |
1237 |
|
1195622 |
zalpha_tr = a_tr_sca(3) |
1238 |
|
|
ELSE |
1239 |
|
1208812 |
zalpha_tr = a_tr_sca(4) |
1240 |
|
|
ENDIF |
1241 |
|
|
|
1242 |
|
2404434 |
zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i)) |
1243 |
|
2404434 |
frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi |
1244 |
|
|
|
1245 |
|
|
! Nucleation with a factor of -1 instead of -0.5 |
1246 |
|
|
zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i)) |
1247 |
|
|
|
1248 |
|
|
ENDIF |
1249 |
|
|
|
1250 |
|
|
ENDDO |
1251 |
|
|
|
1252 |
|
|
! Scavenging through impaction in the underlying layer |
1253 |
|
|
|
1254 |
✓✓ |
224640 |
DO kk = k-1, 1, -1 |
1255 |
|
|
|
1256 |
✓✓ |
212352192 |
DO i = 1, klon |
1257 |
|
|
|
1258 |
✓✓✓✓
|
212340960 |
IF (rneb(i,k).GT.0.0.AND.zprec_cond(i).GT.0.) THEN |
1259 |
|
|
|
1260 |
✓✓ |
24384891 |
IF (t(i,kk) .GE. t_glace_min) THEN |
1261 |
|
19880817 |
zalpha_tr = a_tr_sca(1) |
1262 |
|
|
ELSE |
1263 |
|
4504074 |
zalpha_tr = a_tr_sca(2) |
1264 |
|
|
ENDIF |
1265 |
|
|
|
1266 |
|
24384891 |
zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i)) |
1267 |
|
24384891 |
frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi |
1268 |
|
|
|
1269 |
|
|
ENDIF |
1270 |
|
|
|
1271 |
|
|
ENDDO |
1272 |
|
|
|
1273 |
|
|
ENDDO |
1274 |
|
|
|
1275 |
|
|
!--save some variables for ice supersaturation |
1276 |
|
|
! |
1277 |
✓✓ |
11175840 |
DO i = 1, klon |
1278 |
|
|
! for memory |
1279 |
|
11164608 |
rneb_seri(i,k) = rneb(i,k) |
1280 |
|
|
|
1281 |
|
|
! for diagnostics |
1282 |
|
11164608 |
rnebclr(i,k) = 1.0 - rneb(i,k) - rnebss(i,k) |
1283 |
|
|
|
1284 |
|
11164608 |
qvc(i,k) = zqs(i) * rneb(i,k) |
1285 |
|
11164608 |
qclr(i,k) = MAX(1.e-10,zq(i) - qvc(i,k) - qss(i,k)) !--added by OB because of pathologic cases with the lognormal |
1286 |
|
11175840 |
qcld(i,k) = qvc(i,k) + zcond(i) |
1287 |
|
|
ENDDO |
1288 |
|
|
!q_sat |
1289 |
|
11232 |
CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,1,.false.,qsatl(:,k),zdqs(:)) |
1290 |
|
11520 |
CALL calc_qsat_ecmwf(klon,Tbef(:),qzero(:),pplay(:,k),RTT,2,.false.,qsats(:,k),zdqs(:)) |
1291 |
|
|
|
1292 |
|
|
ENDDO |
1293 |
|
|
|
1294 |
|
|
!====================================================================== |
1295 |
|
|
! END OF VERTICAL LOOP |
1296 |
|
|
!====================================================================== |
1297 |
|
|
|
1298 |
|
|
! Rain or snow at the surface (depending on the first layer temperature) |
1299 |
✓✓ |
286560 |
DO i = 1, klon |
1300 |
|
286272 |
snow(i) = zifl(i) |
1301 |
|
286560 |
rain(i) = zrfl(i) |
1302 |
|
|
! c_iso final output |
1303 |
|
|
ENDDO |
1304 |
|
|
|
1305 |
✗✓ |
288 |
IF (ncoreczq>0) THEN |
1306 |
|
|
WRITE(lunout,*)'WARNING : ZQ in LSCP ',ncoreczq,' val < 1.e-15.' |
1307 |
|
|
ENDIF |
1308 |
|
|
|
1309 |
|
288 |
END SUBROUTINE lscp |
1310 |
|
|
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
1311 |
|
|
|
1312 |
|
|
END MODULE lscp_mod |