1 |
|
|
! $Id: fisrtilp.F90 4535 2023-05-13 12:12:05Z evignon $ |
2 |
|
|
! |
3 |
|
|
! |
4 |
|
|
SUBROUTINE fisrtilp(dtime,paprs,pplay,t,q,ptconv,ratqs, & |
5 |
|
|
d_t, d_q, d_ql, d_qi, rneb, radliq, rain, snow, & |
6 |
|
|
pfrac_impa, pfrac_nucl, pfrac_1nucl, & |
7 |
|
|
frac_impa, frac_nucl, beta, & |
8 |
|
|
prfl, psfl, rhcl, zqta, fraca, & |
9 |
|
|
ztv, zpspsk, ztla, zthl, iflag_cld_th, & |
10 |
|
|
iflag_ice_thermo) |
11 |
|
|
|
12 |
|
|
! |
13 |
|
|
USE dimphy |
14 |
|
|
USE icefrac_lsc_mod ! compute ice fraction (JBM 3/14) |
15 |
|
|
USE print_control_mod, ONLY: prt_level, lunout |
16 |
|
|
USE cloudth_mod |
17 |
|
|
USE ioipsl_getin_p_mod, ONLY : getin_p |
18 |
|
|
USE phys_local_var_mod, ONLY: ql_seri,qs_seri |
19 |
|
|
USE phys_local_var_mod, ONLY: rneblsvol |
20 |
|
|
! flag to include modifications to ensure energy conservation (if flag >0) |
21 |
|
|
USE add_phys_tend_mod, only : fl_cor_ebil |
22 |
|
|
USE lscp_ini_mod, ONLY: iflag_t_glace,t_glace_min, t_glace_max, exposant_glace |
23 |
|
|
USE lscp_ini_mod, ONLY: iflag_cloudth_vert, iflag_rain_incloud_vol |
24 |
|
|
USE lscp_ini_mod, ONLY: coef_eva, coef_eva_i, ffallv_lsc, ffallv_con |
25 |
|
|
USE lscp_ini_mod, ONLY: cld_tau_lsc, cld_tau_con, cld_lc_lsc, cld_lc_con |
26 |
|
|
USE lscp_ini_mod, ONLY: reevap_ice, iflag_bergeron, iflag_fisrtilp_qsat, iflag_pdf |
27 |
|
|
|
28 |
|
|
|
29 |
|
|
IMPLICIT none |
30 |
|
|
!====================================================================== |
31 |
|
|
! Auteur(s): Z.X. Li (LMD/CNRS) |
32 |
|
|
! Date: le 20 mars 1995 |
33 |
|
|
! Objet: condensation et precipitation stratiforme. |
34 |
|
|
! schema de nuage |
35 |
|
|
! Fusion de fisrt (physique sursaturation, P. LeVan K. Laval) |
36 |
|
|
! et ilp (il pleut, L. Li) |
37 |
|
|
! Principales parties: |
38 |
|
|
! P0> Thermalisation des precipitations venant de la couche du dessus |
39 |
|
|
! P1> Evaporation de la precipitation (qui vient du niveau k+1) |
40 |
|
|
! P2> Formation du nuage (en k) |
41 |
|
|
! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau |
42 |
|
|
! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour |
43 |
|
|
! les valeurs de T et Q initiales |
44 |
|
|
! P2.A.2> Prise en compte du couplage entre eau condensee et T. |
45 |
|
|
! P2.A.3> Calcul des valeures finales associees a la formation des nuages |
46 |
|
|
! P2.B> Nuage "tout ou rien" |
47 |
|
|
! P2.C> Prise en compte de la Chaleur latente apres formation nuage |
48 |
|
|
! P3> Formation de la precipitation (en k) |
49 |
|
|
!====================================================================== |
50 |
|
|
! JLD: |
51 |
|
|
! * Routine probablement fausse (au moins incoherente) si thermcep = .false. |
52 |
|
|
! * fl_cor_ebil doit etre > 0 ; |
53 |
|
|
! fl_cor_ebil= 0 pour reproduire anciens bugs |
54 |
|
|
!====================================================================== |
55 |
|
|
include "YOMCST.h" |
56 |
|
|
! |
57 |
|
|
! Principaux inputs: |
58 |
|
|
! |
59 |
|
|
REAL, INTENT(IN) :: dtime ! intervalle du temps (s) |
60 |
|
|
REAL, DIMENSION(klon,klev+1), INTENT(IN) :: paprs ! pression a inter-couche |
61 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: pplay ! pression au milieu de couche |
62 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: t ! temperature (K) |
63 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: q ! humidite specifique (kg/kg) |
64 |
|
|
LOGICAL, DIMENSION(klon,klev), INTENT(IN) :: ptconv ! points ou le schema de conv. prof. est actif |
65 |
|
|
INTEGER, INTENT(IN) :: iflag_cld_th |
66 |
|
|
INTEGER, INTENT(IN) :: iflag_ice_thermo |
67 |
|
|
! |
68 |
|
|
! Inputs lies aux thermiques |
69 |
|
|
! |
70 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: ztv |
71 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: zqta, fraca |
72 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: zpspsk, ztla |
73 |
|
|
REAL, DIMENSION(klon,klev), INTENT(IN) :: zthl |
74 |
|
|
! |
75 |
|
|
! Input/output |
76 |
|
|
REAL, DIMENSION(klon,klev), INTENT(INOUT):: ratqs ! determine la largeur de distribution de vapeur |
77 |
|
|
! |
78 |
|
|
! Principaux outputs: |
79 |
|
|
! |
80 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_t ! incrementation de la temperature (K) |
81 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_q ! incrementation de la vapeur d'eau |
82 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_ql ! incrementation de l'eau liquide |
83 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: d_qi ! incrementation de l'eau glace |
84 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: rneb ! fraction nuageuse |
85 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: radliq ! eau liquide utilisee dans rayonnements |
86 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: rhcl ! humidite relative en ciel clair |
87 |
|
|
REAL, DIMENSION(klon), INTENT(OUT) :: rain |
88 |
|
|
REAL, DIMENSION(klon), INTENT(OUT) :: snow |
89 |
|
|
REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: prfl |
90 |
|
|
REAL, DIMENSION(klon,klev+1), INTENT(OUT) :: psfl |
91 |
|
|
|
92 |
|
|
!AA |
93 |
|
|
! Coeffients de fraction lessivee : pour OFF-LINE |
94 |
|
|
! |
95 |
|
|
REAL, DIMENSION(klon,klev), INTENT(INOUT) :: pfrac_nucl |
96 |
|
|
REAL, DIMENSION(klon,klev), INTENT(INOUT) :: pfrac_1nucl |
97 |
|
|
REAL, DIMENSION(klon,klev), INTENT(INOUT) :: pfrac_impa |
98 |
|
|
! |
99 |
|
|
! Fraction d'aerosols lessivee par impaction et par nucleation |
100 |
|
|
! POur ON-LINE |
101 |
|
|
! |
102 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_impa |
103 |
|
|
REAL, DIMENSION(klon,klev), INTENT(OUT) :: frac_nucl |
104 |
|
|
!AA |
105 |
|
|
! -------------------------------------------------------------------------------- |
106 |
|
|
! |
107 |
|
|
! Options du programme: |
108 |
|
|
! |
109 |
|
|
REAL, SAVE :: seuil_neb=0.001 ! un nuage existe vraiment au-dela |
110 |
|
|
!$OMP THREADPRIVATE(seuil_neb) |
111 |
|
|
|
112 |
|
|
!<LTP |
113 |
|
|
REAL smallestreal |
114 |
|
|
REAL, SAVE :: rain_int_min=0.001 !intensité locale minimum pour la pluie avant diminution de la fraction précipitante associée = 0.001 mm/s |
115 |
|
|
!>LTP |
116 |
|
|
!$OMP THREADPRIVATE(rain_int_min) |
117 |
|
|
|
118 |
|
|
|
119 |
|
|
INTEGER ninter ! sous-intervals pour la precipitation |
120 |
|
|
PARAMETER (ninter=5) |
121 |
|
|
INTEGER,SAVE :: iflag_evap_prec=1 ! evaporation de la pluie |
122 |
|
|
!$OMP THREADPRIVATE(iflag_evap_prec) |
123 |
|
|
! |
124 |
|
|
LOGICAL cpartiel ! condensation partielle |
125 |
|
|
PARAMETER (cpartiel=.TRUE.) |
126 |
|
|
REAL t_coup |
127 |
|
|
PARAMETER (t_coup=234.0) |
128 |
|
|
REAL DDT0 |
129 |
|
|
PARAMETER (DDT0=.01) |
130 |
|
|
REAL ztfondue |
131 |
|
|
PARAMETER (ztfondue=278.15) |
132 |
|
|
! -------------------------------------------------------------------------------- |
133 |
|
|
! |
134 |
|
|
! Variables locales: |
135 |
|
|
! |
136 |
|
|
INTEGER i, k, n, kk |
137 |
|
|
INTEGER,save::itap=0 |
138 |
|
|
!$OMP THREADPRIVATE(itap) |
139 |
|
|
|
140 |
|
|
REAL qsl, qsi |
141 |
|
|
real zct ,zcl |
142 |
|
|
INTEGER ncoreczq |
143 |
|
|
REAL ctot(klon,klev) |
144 |
|
|
REAL ctot_vol(klon,klev) |
145 |
|
|
REAL zqs(klon), zdqs(klon), zdelta, zcor, zcvm5 |
146 |
|
|
REAL zdqsdT_raw(klon) |
147 |
|
|
REAL Tbef(klon),qlbef(klon),DT(klon),num,denom |
148 |
|
|
|
149 |
|
|
logical lognormale(klon) |
150 |
|
|
logical ice_thermo |
151 |
|
|
LOGICAL convergence(klon) |
152 |
|
|
INTEGER n_i(klon), iter |
153 |
|
|
REAL cste |
154 |
|
|
|
155 |
|
|
real zpdf_sig(klon),zpdf_k(klon),zpdf_delta(klon) |
156 |
|
|
real Zpdf_a(klon),zpdf_b(klon),zpdf_e1(klon),zpdf_e2(klon) |
157 |
|
|
real erf |
158 |
|
|
REAL qcloud(klon) |
159 |
|
|
|
160 |
|
|
REAL zrfl(klon), zrfln(klon), zqev, zqevt |
161 |
|
|
!<LTP |
162 |
|
|
REAL zrflclr(klon), zrflcld(klon) |
163 |
|
|
REAL d_zrfl_clr_cld(klon), d_zifl_clr_cld(klon) |
164 |
|
|
REAL d_zrfl_cld_clr(klon), d_zifl_cld_clr(klon) |
165 |
|
|
!>LTP |
166 |
|
|
|
167 |
|
|
REAL zifl(klon), zifln(klon), zqev0,zqevi, zqevti |
168 |
|
|
!<LTP |
169 |
|
|
REAL ziflclr(klon), ziflcld(klon) |
170 |
|
|
!>LTP |
171 |
|
|
REAL zoliq(klon), zcond(klon), zq(klon), zqn(klon), zdelq |
172 |
|
|
REAL zoliqp(klon), zoliqi(klon) |
173 |
|
|
REAL zt(klon) |
174 |
|
|
! JBM (3/14) nexpo is replaced by exposant_glace |
175 |
|
|
! REAL nexpo ! exponentiel pour glace/eau |
176 |
|
|
! INTEGER, PARAMETER :: nexpo=6 |
177 |
|
|
INTEGER exposant_glace_old |
178 |
|
|
REAL t_glace_min_old |
179 |
|
|
REAL zdz(klon),zrho(klon),ztot , zrhol(klon) |
180 |
|
|
REAL zchau ,zfroi ,zfice(klon),zneb(klon),znebprecip(klon) |
181 |
|
|
!<LTP |
182 |
|
|
REAL znebprecipclr(klon), znebprecipcld(klon) |
183 |
|
|
REAL tot_zneb(klon), tot_znebn(klon), d_tot_zneb(klon) |
184 |
|
|
REAL d_znebprecip_clr_cld(klon), d_znebprecip_cld_clr(klon) |
185 |
|
|
!>LTP |
186 |
|
|
|
187 |
|
|
REAL zmelt, zpluie, zice |
188 |
|
|
REAL dzfice(klon) |
189 |
|
|
REAL zsolid |
190 |
|
|
!!!! |
191 |
|
|
! Variables pour Bergeron |
192 |
|
|
REAL zcp, coef1, DeltaT, Deltaq, Deltaqprecl |
193 |
|
|
REAL zqpreci(klon), zqprecl(klon) |
194 |
|
|
! Variable pour conservation enegie des precipitations |
195 |
|
|
REAL zmqc(klon) |
196 |
|
|
! |
197 |
|
|
LOGICAL appel1er |
198 |
|
|
SAVE appel1er |
199 |
|
|
!$OMP THREADPRIVATE(appel1er) |
200 |
|
|
! |
201 |
|
|
! iflag_oldbug_fisrtilp=0 enleve le BUG par JYG : tglace_min -> tglace_max |
202 |
|
|
! iflag_oldbug_fisrtilp=1 ajoute le BUG |
203 |
|
|
INTEGER,SAVE :: iflag_oldbug_fisrtilp=0 !=0 sans bug |
204 |
|
|
!$OMP THREADPRIVATE(iflag_oldbug_fisrtilp) |
205 |
|
|
!--------------------------------------------------------------- |
206 |
|
|
! |
207 |
|
|
!AA Variables traceurs: |
208 |
|
|
!AA Provisoire !!! Parametres alpha du lessivage |
209 |
|
|
!AA A priori on a 4 scavenging # possibles |
210 |
|
|
! |
211 |
|
|
REAL a_tr_sca(4) |
212 |
|
|
save a_tr_sca |
213 |
|
|
!$OMP THREADPRIVATE(a_tr_sca) |
214 |
|
|
! |
215 |
|
|
! Variables intermediaires |
216 |
|
|
! |
217 |
|
|
REAL zalpha_tr |
218 |
|
|
REAL zfrac_lessi |
219 |
|
|
REAL zprec_cond(klon) |
220 |
|
|
!AA |
221 |
|
|
! RomP >>> 15 nov 2012 |
222 |
|
|
REAL beta(klon,klev) ! taux de conversion de l'eau cond |
223 |
|
|
! RomP <<< |
224 |
|
|
REAL zmair(klon), zcpair, zcpeau |
225 |
|
|
! Pour la conversion eau-neige |
226 |
|
|
REAL zlh_solid(klon), zm_solid |
227 |
|
|
!--------------------------------------------------------------- |
228 |
|
|
! |
229 |
|
|
! Fonctions en ligne: |
230 |
|
|
! |
231 |
|
|
REAL fallvs,fallvc ! Vitesse de chute pour cristaux de glace |
232 |
|
|
! (Heymsfield & Donner, 1990) |
233 |
|
|
REAL zzz |
234 |
|
|
include "YOETHF.h" |
235 |
|
|
include "FCTTRE.h" |
236 |
|
|
fallvc (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_con |
237 |
|
|
fallvs (zzz) = 3.29/2.0 * ((zzz)**0.16) * ffallv_lsc |
238 |
|
|
! |
239 |
|
|
DATA appel1er /.TRUE./ |
240 |
|
|
!ym |
241 |
|
|
!CR: pour iflag_ice_thermo=2, on active que la convection |
242 |
|
|
! ice_thermo = iflag_ice_thermo .GE. 1 |
243 |
|
|
|
244 |
|
|
|
245 |
|
|
itap=itap+1 |
246 |
|
|
znebprecip(:)=0. |
247 |
|
|
|
248 |
|
|
!<LTP |
249 |
|
|
smallestreal=1.e-9 |
250 |
|
|
znebprecipclr(:)=0. |
251 |
|
|
znebprecipcld(:)=0. |
252 |
|
|
!>LTP |
253 |
|
|
|
254 |
|
|
ice_thermo = (iflag_ice_thermo .EQ. 1).OR.(iflag_ice_thermo .GE. 3) |
255 |
|
|
zdelq=0.0 |
256 |
|
|
ctot_vol(1:klon,1:klev)=0.0 |
257 |
|
|
rneblsvol(1:klon,1:klev)=0.0 |
258 |
|
|
|
259 |
|
|
if (prt_level>9)write(lunout,*)'NUAGES4 A. JAM' |
260 |
|
|
IF (appel1er) THEN |
261 |
|
|
CALL getin_p('iflag_oldbug_fisrtilp',iflag_oldbug_fisrtilp) |
262 |
|
|
CALL getin_p('iflag_evap_prec',iflag_evap_prec) |
263 |
|
|
CALL getin_p('seuil_neb',seuil_neb) |
264 |
|
|
!<LTP |
265 |
|
|
CALL getin_p('rain_int_min',rain_int_min) |
266 |
|
|
!>LTP |
267 |
|
|
write(lunout,*)' iflag_oldbug_fisrtilp =',iflag_oldbug_fisrtilp |
268 |
|
|
! |
269 |
|
|
WRITE(lunout,*) 'fisrtilp, ninter:', ninter |
270 |
|
|
WRITE(lunout,*) 'fisrtilp, iflag_evap_prec:', iflag_evap_prec |
271 |
|
|
!<LTP |
272 |
|
|
WRITE(lunout,*) 'fisrtilp, rain_int_min:', rain_int_min |
273 |
|
|
!>LTP |
274 |
|
|
WRITE(lunout,*) 'fisrtilp, cpartiel:', cpartiel |
275 |
|
|
WRITE(lunout,*) 'FISRTILP VERSION LUDO' |
276 |
|
|
|
277 |
|
|
IF (ABS(dtime/REAL(ninter)-360.0).GT.0.001) THEN |
278 |
|
|
WRITE(lunout,*) 'fisrtilp: Ce n est pas prevu, voir Z.X.Li', dtime |
279 |
|
|
WRITE(lunout,*) 'Je prefere un sous-intervalle de 6 minutes' |
280 |
|
|
! CALL abort |
281 |
|
|
ENDIF |
282 |
|
|
appel1er = .FALSE. |
283 |
|
|
! |
284 |
|
|
!AA initialiation provisoire |
285 |
|
|
a_tr_sca(1) = -0.5 |
286 |
|
|
a_tr_sca(2) = -0.5 |
287 |
|
|
a_tr_sca(3) = -0.5 |
288 |
|
|
a_tr_sca(4) = -0.5 |
289 |
|
|
! |
290 |
|
|
!AA Initialisation a 1 des coefs des fractions lessivees |
291 |
|
|
! |
292 |
|
|
!cdir collapse |
293 |
|
|
DO k = 1, klev |
294 |
|
|
DO i = 1, klon |
295 |
|
|
pfrac_nucl(i,k)=1. |
296 |
|
|
pfrac_1nucl(i,k)=1. |
297 |
|
|
pfrac_impa(i,k)=1. |
298 |
|
|
beta(i,k)=0. !RomP initialisation |
299 |
|
|
ENDDO |
300 |
|
|
ENDDO |
301 |
|
|
|
302 |
|
|
ENDIF ! test sur appel1er |
303 |
|
|
! |
304 |
|
|
!MAf Initialisation a 0 de zoliq |
305 |
|
|
! DO i = 1, klon |
306 |
|
|
! zoliq(i)=0. |
307 |
|
|
! ENDDO |
308 |
|
|
! Determiner les nuages froids par leur temperature |
309 |
|
|
! nexpo regle la raideur de la transition eau liquide / eau glace. |
310 |
|
|
! |
311 |
|
|
!CR: on est oblige de definir des valeurs fisrt car les valeurs de newmicro ne sont pas les memes par defaut |
312 |
|
|
IF (iflag_t_glace.EQ.0) THEN |
313 |
|
|
! ztglace = RTT - 15.0 |
314 |
|
|
t_glace_min_old = RTT - 15.0 |
315 |
|
|
!AJ< |
316 |
|
|
IF (ice_thermo) THEN |
317 |
|
|
! nexpo = 2 |
318 |
|
|
exposant_glace_old = 2 |
319 |
|
|
ELSE |
320 |
|
|
! nexpo = 6 |
321 |
|
|
exposant_glace_old = 6 |
322 |
|
|
ENDIF |
323 |
|
|
|
324 |
|
|
ENDIF |
325 |
|
|
|
326 |
|
|
!! RLVTT = 2.501e6 ! pas de redefinition des constantes physiques (jyg) |
327 |
|
|
!! RLSTT = 2.834e6 ! pas de redefinition des constantes physiques (jyg) |
328 |
|
|
!>AJ |
329 |
|
|
!cc nexpo = 1 |
330 |
|
|
! |
331 |
|
|
! Initialiser les sorties: |
332 |
|
|
! |
333 |
|
|
!cdir collapse |
334 |
|
|
DO k = 1, klev+1 |
335 |
|
|
DO i = 1, klon |
336 |
|
|
prfl(i,k) = 0.0 |
337 |
|
|
psfl(i,k) = 0.0 |
338 |
|
|
ENDDO |
339 |
|
|
ENDDO |
340 |
|
|
|
341 |
|
|
!cdir collapse |
342 |
|
|
|
343 |
|
|
DO k = 1, klev |
344 |
|
|
DO i = 1, klon |
345 |
|
|
d_t(i,k) = 0.0 |
346 |
|
|
d_q(i,k) = 0.0 |
347 |
|
|
d_ql(i,k) = 0.0 |
348 |
|
|
d_qi(i,k) = 0.0 |
349 |
|
|
rneb(i,k) = 0.0 |
350 |
|
|
radliq(i,k) = 0.0 |
351 |
|
|
frac_nucl(i,k) = 1. |
352 |
|
|
frac_impa(i,k) = 1. |
353 |
|
|
ENDDO |
354 |
|
|
ENDDO |
355 |
|
|
DO i = 1, klon |
356 |
|
|
rain(i) = 0.0 |
357 |
|
|
snow(i) = 0.0 |
358 |
|
|
zoliq(i)=0. |
359 |
|
|
! ENDDO |
360 |
|
|
! |
361 |
|
|
! Initialiser le flux de precipitation a zero |
362 |
|
|
! |
363 |
|
|
! DO i = 1, klon |
364 |
|
|
zrfl(i) = 0.0 |
365 |
|
|
zifl(i) = 0.0 |
366 |
|
|
!<LTP |
367 |
|
|
zrflclr(i) = 0.0 |
368 |
|
|
ziflclr(i) = 0.0 |
369 |
|
|
zrflcld(i) = 0.0 |
370 |
|
|
ziflcld(i) = 0.0 |
371 |
|
|
tot_zneb(i) = 0.0 |
372 |
|
|
tot_znebn(i) = 0.0 |
373 |
|
|
d_tot_zneb(i) = 0.0 |
374 |
|
|
!>LTP |
375 |
|
|
|
376 |
|
|
zneb(i) = seuil_neb |
377 |
|
|
ENDDO |
378 |
|
|
! |
379 |
|
|
! |
380 |
|
|
!AA Pour plus de securite |
381 |
|
|
|
382 |
|
|
zalpha_tr = 0. |
383 |
|
|
zfrac_lessi = 0. |
384 |
|
|
|
385 |
|
|
!AA================================================================== |
386 |
|
|
! |
387 |
|
|
ncoreczq=0 |
388 |
|
|
! BOUCLE VERTICALE (DU HAUT VERS LE BAS) |
389 |
|
|
! |
390 |
|
|
DO k = klev, 1, -1 |
391 |
|
|
! |
392 |
|
|
!AA=============================================================== |
393 |
|
|
! |
394 |
|
|
! Initialisation temperature et vapeur |
395 |
|
|
DO i = 1, klon |
396 |
|
|
zt(i)=t(i,k) |
397 |
|
|
zq(i)=q(i,k) |
398 |
|
|
ENDDO |
399 |
|
|
! |
400 |
|
|
! ---------------------------------------------------------------- |
401 |
|
|
! P0> Thermalisation des precipitations venant de la couche du dessus |
402 |
|
|
! ---------------------------------------------------------------- |
403 |
|
|
! Calculer la varition de temp. de l'air du a la chaleur sensible |
404 |
|
|
! transporter par la pluie. On thermalise la pluie avec l'air de la couche. |
405 |
|
|
! Cette quantite de pluie qui est thermalisee, et devra continue a l'etre lors |
406 |
|
|
! des differentes transformations thermodynamiques. Cette masse d'eau doit |
407 |
|
|
! donc etre ajoute a l'humidite de la couche lorsque l'on calcule la variation |
408 |
|
|
! de l'enthalpie de la couche avec la temperature |
409 |
|
|
! Variables calculees ou modifiees: |
410 |
|
|
! - zt: temperature de la cocuhe |
411 |
|
|
! - zmqc: masse de precip qui doit etre thermalisee |
412 |
|
|
! |
413 |
|
|
IF(k.LE.klevm1) THEN |
414 |
|
|
DO i = 1, klon |
415 |
|
|
!IM |
416 |
|
|
zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG |
417 |
|
|
! il n'y a pas encore d'eau liquide ni glace dans la maiille, donc zq suffit |
418 |
|
|
zcpair=RCPD*(1.0+RVTMP2*zq(i)) |
419 |
|
|
zcpeau=RCPD*RVTMP2 |
420 |
|
|
if (fl_cor_ebil .GT. 0) then |
421 |
|
|
! zmqc: masse de precip qui doit etre thermalisee avec l'air de la couche atm |
422 |
|
|
! pour s'assurer que la precip arrivant au sol aura bien la temperature de la |
423 |
|
|
! derniere couche |
424 |
|
|
zmqc(i) = (zrfl(i)+zifl(i))*dtime/zmair(i) |
425 |
|
|
! t(i,k+1)+d_t(i,k+1): nouvelle temp de la couche au dessus |
426 |
|
|
zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zmqc(i)*zcpeau + zcpair*zt(i) ) & |
427 |
|
|
/ (zcpair + zmqc(i)*zcpeau) |
428 |
|
|
else ! si on maintient les anciennes erreurs |
429 |
|
|
zt(i) = ( (t(i,k+1)+d_t(i,k+1))*zrfl(i)*dtime*zcpeau & |
430 |
|
|
+ zmair(i)*zcpair*zt(i) ) & |
431 |
|
|
/ (zmair(i)*zcpair + zrfl(i)*dtime*zcpeau) |
432 |
|
|
end if |
433 |
|
|
ENDDO |
434 |
|
|
ELSE ! IF(k.LE.klevm1) |
435 |
|
|
DO i = 1, klon |
436 |
|
|
zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG |
437 |
|
|
zmqc(i) = 0. |
438 |
|
|
ENDDO |
439 |
|
|
ENDIF ! end IF(k.LE.klevm1) |
440 |
|
|
! |
441 |
|
|
! ---------------------------------------------------------------- |
442 |
|
|
! P1> Calcul de l'evaporation de la precipitation |
443 |
|
|
! ---------------------------------------------------------------- |
444 |
|
|
! On evapore une partie des precipitations venant de la maille du dessus. |
445 |
|
|
! On calcule l'evaporation et la sublimation des precipitations, jusqu'a |
446 |
|
|
! ce que la fraction de cette couche qui est sous le nuage soit saturee. |
447 |
|
|
! Variables calculees ou modifiees: |
448 |
|
|
! - zrfl et zifl: flux de precip liquide et glace |
449 |
|
|
! - zq, zt: humidite et temperature de la cocuhe |
450 |
|
|
! - zmqc: masse de precip qui doit etre thermalisee |
451 |
|
|
! |
452 |
|
|
IF (iflag_evap_prec>=1) THEN |
453 |
|
|
DO i = 1, klon |
454 |
|
|
! S'il y a des precipitations |
455 |
|
|
IF (zrfl(i)+zifl(i).GT.0.) THEN |
456 |
|
|
! Calcul du qsat |
457 |
|
|
IF (thermcep) THEN |
458 |
|
|
zdelta=MAX(0.,SIGN(1.,RTT-zt(i))) |
459 |
|
|
zqs(i)= R2ES*FOEEW(zt(i),zdelta)/pplay(i,k) |
460 |
|
|
zqs(i)=MIN(0.5,zqs(i)) |
461 |
|
|
zcor=1./(1.-RETV*zqs(i)) |
462 |
|
|
zqs(i)=zqs(i)*zcor |
463 |
|
|
ELSE |
464 |
|
|
IF (zt(i) .LT. t_coup) THEN |
465 |
|
|
zqs(i) = qsats(zt(i)) / pplay(i,k) |
466 |
|
|
ELSE |
467 |
|
|
zqs(i) = qsatl(zt(i)) / pplay(i,k) |
468 |
|
|
ENDIF |
469 |
|
|
ENDIF |
470 |
|
|
ENDIF ! (zrfl(i)+zifl(i).GT.0.) |
471 |
|
|
ENDDO |
472 |
|
|
!AJ< |
473 |
|
|
|
474 |
|
|
IF (.NOT. ice_thermo) THEN |
475 |
|
|
DO i = 1, klon |
476 |
|
|
! S'il y a des precipitations |
477 |
|
|
IF (zrfl(i)+zifl(i).GT.0.) THEN |
478 |
|
|
! Evap max pour ne pas saturer la fraction sous le nuage |
479 |
|
|
! Evap max jusqu'à atteindre la saturation dans la partie |
480 |
|
|
! de la maille qui est sous le nuage de la couche du dessus |
481 |
|
|
!!! On ne tient compte de cette fraction que sous une seule |
482 |
|
|
!!! couche sous le nuage |
483 |
|
|
zqev = MAX (0.0, (zqs(i)-zq(i))*zneb(i) ) |
484 |
|
|
! Ajout de la prise en compte des precip a thermiser |
485 |
|
|
! avec petite reecriture |
486 |
|
|
if (fl_cor_ebil .GT. 0) then ! nouveau |
487 |
|
|
! Calcul de l'evaporation du flux de precip herite |
488 |
|
|
! d'au-dessus |
489 |
|
|
zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) & |
490 |
|
|
* zmair(i)/pplay(i,k)*zt(i)*RD |
491 |
|
|
zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) * dtime/zmair(i) |
492 |
|
|
|
493 |
|
|
! Seuil pour ne pas saturer la fraction sous le nuage |
494 |
|
|
zqev = MIN (zqev, zqevt) |
495 |
|
|
! Nouveau flux de precip |
496 |
|
|
zrfln(i) = zrfl(i) - zqev*zmair(i)/dtime |
497 |
|
|
! Aucun flux liquide pour T < t_coup, on reevapore tout. |
498 |
|
|
IF (zt(i) .LT. t_coup.and.reevap_ice) THEN |
499 |
|
|
zrfln(i)=0. |
500 |
|
|
zqev = (zrfl(i)-zrfln(i))/zmair(i)*dtime |
501 |
|
|
END IF |
502 |
|
|
! Nouvelle vapeur |
503 |
|
|
zq(i) = zq(i) + zqev |
504 |
|
|
zmqc(i) = zmqc(i)-zqev |
505 |
|
|
! Nouvelle temperature (chaleur latente) |
506 |
|
|
zt(i) = zt(i) - zqev & |
507 |
|
|
* RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) |
508 |
|
|
!!JLD debut de partie a supprimer a terme |
509 |
|
|
else ! if (fl_cor_ebil .GT. 0) |
510 |
|
|
! Calcul de l'evaporation du flux de precip herite |
511 |
|
|
! d'au-dessus |
512 |
|
|
zqevt = coef_eva * (1.0-zq(i)/zqs(i)) * SQRT(zrfl(i)) & |
513 |
|
|
* (paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG |
514 |
|
|
zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) & |
515 |
|
|
* RG*dtime/(paprs(i,k)-paprs(i,k+1)) |
516 |
|
|
! Seuil pour ne pas saturer la fraction sous le nuage |
517 |
|
|
zqev = MIN (zqev, zqevt) |
518 |
|
|
! Nouveau flux de precip |
519 |
|
|
zrfln(i) = zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) & |
520 |
|
|
/RG/dtime |
521 |
|
|
! Aucun flux liquide pour T < t_coup |
522 |
|
|
IF (zt(i) .LT. t_coup.and.reevap_ice) zrfln(i)=0. |
523 |
|
|
! Nouvelle vapeur |
524 |
|
|
zq(i) = zq(i) - (zrfln(i)-zrfl(i)) & |
525 |
|
|
* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime |
526 |
|
|
! Nouvelle temperature (chaleur latente) |
527 |
|
|
zt(i) = zt(i) + (zrfln(i)-zrfl(i)) & |
528 |
|
|
* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime & |
529 |
|
|
* RLVTT/RCPD/(1.0+RVTMP2*zq(i)) |
530 |
|
|
end if ! if (fl_cor_ebil .GT. 0) |
531 |
|
|
!!JLD fin de partie a supprimer a terme |
532 |
|
|
zrfl(i) = zrfln(i) |
533 |
|
|
zifl(i) = 0. |
534 |
|
|
ENDIF ! (zrfl(i)+zifl(i).GT.0.) |
535 |
|
|
ENDDO |
536 |
|
|
! |
537 |
|
|
ELSE ! (.NOT. ice_thermo) |
538 |
|
|
! ================================ |
539 |
|
|
! Avec thermodynamique de la glace |
540 |
|
|
! ================================ |
541 |
|
|
DO i = 1, klon |
542 |
|
|
|
543 |
|
|
|
544 |
|
|
!AJ< |
545 |
|
|
! S'il y a des precipitations |
546 |
|
|
IF (zrfl(i)+zifl(i).GT.0.) THEN |
547 |
|
|
|
548 |
|
|
!LTP< |
549 |
|
|
!On ne tient compte que du flux de précipitation en ciel clair dans le calcul de l'évaporation. |
550 |
|
|
IF (iflag_evap_prec==4) THEN |
551 |
|
|
zrfl(i) = zrflclr(i) |
552 |
|
|
zifl(i) = ziflclr(i) |
553 |
|
|
ENDIF |
554 |
|
|
|
555 |
|
|
!>LTP |
556 |
|
|
|
557 |
|
|
IF (iflag_evap_prec==1) THEN |
558 |
|
|
znebprecip(i)=zneb(i) |
559 |
|
|
ELSE |
560 |
|
|
znebprecip(i)=MAX(zneb(i),znebprecip(i)) |
561 |
|
|
ENDIF |
562 |
|
|
|
563 |
|
|
IF (iflag_evap_prec==4) THEN |
564 |
|
|
! Evap max pour ne pas saturer toute la maille |
565 |
|
|
zqev0 = MAX (0.0, zqs(i)-zq(i)) |
566 |
|
|
ELSE |
567 |
|
|
! Evap max pour ne pas saturer la fraction sous le nuage |
568 |
|
|
zqev0 = MAX (0.0, (zqs(i)-zq(i))*znebprecip(i) ) |
569 |
|
|
ENDIF |
570 |
|
|
|
571 |
|
|
!JAM |
572 |
|
|
! On differencie qsat pour l'eau et la glace |
573 |
|
|
! Si zdelta=1. --> glace |
574 |
|
|
! Si zdelta=0. --> eau liquide |
575 |
|
|
|
576 |
|
|
! Calcul du qsat par rapport a l'eau liquide |
577 |
|
|
qsl= R2ES*FOEEW(zt(i),0.)/pplay(i,k) |
578 |
|
|
qsl= MIN(0.5,qsl) |
579 |
|
|
zcor= 1./(1.-RETV*qsl) |
580 |
|
|
qsl= qsl*zcor |
581 |
|
|
|
582 |
|
|
! Calcul de l'evaporation du flux de precip venant du dessus |
583 |
|
|
! Formulation en racine du flux de precip |
584 |
|
|
! (Klemp & Wilhelmson, 1978; Sundqvist, 1988) |
585 |
|
|
IF (iflag_evap_prec==3) THEN |
586 |
|
|
zqevt = znebprecip(i)*coef_eva*(1.0-zq(i)/qsl) & |
587 |
|
|
*SQRT(zrfl(i)/max(1.e-4,znebprecip(i))) & |
588 |
|
|
*(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG |
589 |
|
|
!<LTP |
590 |
|
|
ELSE IF (iflag_evap_prec==4) THEN |
591 |
|
|
zqevt = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsl) & |
592 |
|
|
*SQRT(zrfl(i)/max(1.e-8,znebprecipclr(i))) & |
593 |
|
|
*(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG |
594 |
|
|
!>LTP |
595 |
|
|
ELSE |
596 |
|
|
zqevt = 1.*coef_eva*(1.0-zq(i)/qsl)*SQRT(zrfl(i)) & |
597 |
|
|
*(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG |
598 |
|
|
ENDIF |
599 |
|
|
|
600 |
|
|
|
601 |
|
|
zqevt = MAX(0.0,MIN(zqevt,zrfl(i))) & |
602 |
|
|
*RG*dtime/(paprs(i,k)-paprs(i,k+1)) |
603 |
|
|
|
604 |
|
|
! Calcul du qsat par rapport a la glace |
605 |
|
|
qsi= R2ES*FOEEW(zt(i),1.)/pplay(i,k) |
606 |
|
|
qsi= MIN(0.5,qsi) |
607 |
|
|
zcor= 1./(1.-RETV*qsi) |
608 |
|
|
qsi= qsi*zcor |
609 |
|
|
|
610 |
|
|
! Calcul de la sublimation du flux de precip solide herite |
611 |
|
|
! d'au-dessus |
612 |
|
|
IF (iflag_evap_prec==3) THEN |
613 |
|
|
zqevti = znebprecip(i)*coef_eva*(1.0-zq(i)/qsi) & |
614 |
|
|
*SQRT(zifl(i)/max(1.e-4,znebprecip(i))) & |
615 |
|
|
*(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG |
616 |
|
|
!<LTP |
617 |
|
|
ELSE IF (iflag_evap_prec==4) THEN |
618 |
|
|
zqevti = znebprecipclr(i)*coef_eva*(1.0-zq(i)/qsi) & |
619 |
|
|
*SQRT(zifl(i)/max(1.e-8,znebprecipclr(i))) & |
620 |
|
|
*(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG |
621 |
|
|
!>LTP |
622 |
|
|
ELSE |
623 |
|
|
zqevti = 1.*coef_eva*(1.0-zq(i)/qsi)*SQRT(zifl(i)) & |
624 |
|
|
*(paprs(i,k)-paprs(i,k+1))/pplay(i,k)*zt(i)*RD/RG |
625 |
|
|
ENDIF |
626 |
|
|
zqevti = MAX(0.0,MIN(zqevti,zifl(i))) & |
627 |
|
|
*RG*dtime/(paprs(i,k)-paprs(i,k+1)) |
628 |
|
|
|
629 |
|
|
|
630 |
|
|
!JAM |
631 |
|
|
! Limitation de l'evaporation. On s'assure qu'on ne sature pas |
632 |
|
|
! la fraction de la couche sous le nuage sinon on repartit zqev0 |
633 |
|
|
! en conservant la proportion liquide / glace |
634 |
|
|
|
635 |
|
|
IF (zqevt+zqevti.GT.zqev0) THEN |
636 |
|
|
zqev=zqev0*zqevt/(zqevt+zqevti) |
637 |
|
|
zqevi=zqev0*zqevti/(zqevt+zqevti) |
638 |
|
|
ELSE |
639 |
|
|
!JLD je ne comprends pas les lignes ci-dessous. On repartit les precips |
640 |
|
|
! liquides et solides meme si on ne sature pas la couche. |
641 |
|
|
! A mon avis, le test est inutile, et il faudrait tout remplacer par: |
642 |
|
|
! zqev=zqevt |
643 |
|
|
! zqevi=zqevti |
644 |
|
|
IF (zqevt+zqevti.GT.0.) THEN |
645 |
|
|
zqev=MIN(zqev0*zqevt/(zqevt+zqevti),zqevt) |
646 |
|
|
zqevi=MIN(zqev0*zqevti/(zqevt+zqevti),zqevti) |
647 |
|
|
ELSE |
648 |
|
|
zqev=0. |
649 |
|
|
zqevi=0. |
650 |
|
|
ENDIF |
651 |
|
|
ENDIF |
652 |
|
|
|
653 |
|
|
! Nouveaux flux de precip liquide et solide |
654 |
|
|
zrfln(i) = Max(0.,zrfl(i) - zqev*(paprs(i,k)-paprs(i,k+1)) & |
655 |
|
|
/RG/dtime) |
656 |
|
|
zifln(i) = Max(0.,zifl(i) - zqevi*(paprs(i,k)-paprs(i,k+1)) & |
657 |
|
|
/RG/dtime) |
658 |
|
|
|
659 |
|
|
! Mise a jour de la vapeur, temperature et flux de precip |
660 |
|
|
zq(i) = zq(i) - (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) & |
661 |
|
|
* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime |
662 |
|
|
if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips |
663 |
|
|
zmqc(i) = zmqc(i) + (zrfln(i)+zifln(i)-zrfl(i)-zifl(i)) & |
664 |
|
|
* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime |
665 |
|
|
zt(i) = zt(i) + (zrfln(i)-zrfl(i)) & |
666 |
|
|
* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime & |
667 |
|
|
* RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) & |
668 |
|
|
+ (zifln(i)-zifl(i)) & |
669 |
|
|
* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime & |
670 |
|
|
* RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) |
671 |
|
|
else ! sans correction thermalisation des precips |
672 |
|
|
zt(i) = zt(i) + (zrfln(i)-zrfl(i)) & |
673 |
|
|
* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime & |
674 |
|
|
* RLVTT/RCPD/(1.0+RVTMP2*zq(i)) & |
675 |
|
|
+ (zifln(i)-zifl(i)) & |
676 |
|
|
* (RG/(paprs(i,k)-paprs(i,k+1)))*dtime & |
677 |
|
|
* RLSTT/RCPD/(1.0+RVTMP2*zq(i)) |
678 |
|
|
end if |
679 |
|
|
! Nouvelles vaeleurs des precips liquides et solides |
680 |
|
|
zrfl(i) = zrfln(i) |
681 |
|
|
zifl(i) = zifln(i) |
682 |
|
|
|
683 |
|
|
!<LTP |
684 |
|
|
IF (iflag_evap_prec==4) THEN |
685 |
|
|
zrflclr(i) = zrfl(i) |
686 |
|
|
ziflclr(i) = zifl(i) |
687 |
|
|
IF(zrflclr(i) + ziflclr(i) .LE. 0) THEN |
688 |
|
|
znebprecipclr(i) = 0. |
689 |
|
|
ENDIF |
690 |
|
|
zrfl(i) = zrflclr(i) + zrflcld(i) |
691 |
|
|
zifl(i) = ziflclr(i) + ziflcld(i) |
692 |
|
|
ENDIF |
693 |
|
|
!>LTP |
694 |
|
|
|
695 |
|
|
|
696 |
|
|
! print*,'REEVAP ',itap,k,znebprecip(1),zqev0,zqev,zqevi,zrfl(1) |
697 |
|
|
|
698 |
|
|
!CR ATTENTION: deplacement de la fonte de la glace |
699 |
|
|
!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg |
700 |
|
|
!!! zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2 !!!!!!!!! jyg |
701 |
|
|
!jyg : Bug !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! jyg |
702 |
|
|
zmelt = ((zt(i)-273.15)/(ztfondue-273.15)) ! jyg |
703 |
|
|
! fraction de la precip solide qui est fondue |
704 |
|
|
zmelt = MIN(MAX(zmelt,0.),1.) |
705 |
|
|
! Fusion de la glace |
706 |
|
|
!<LTP |
707 |
|
|
IF (iflag_evap_prec==4) THEN |
708 |
|
|
zrflclr(i)=zrflclr(i)+zmelt*ziflclr(i) |
709 |
|
|
zrflcld(i)=zrflcld(i)+zmelt*ziflcld(i) |
710 |
|
|
zrfl(i)=zrflclr(i)+zrflcld(i) |
711 |
|
|
!>LTP |
712 |
|
|
ELSE |
713 |
|
|
zrfl(i)=zrfl(i)+zmelt*zifl(i) |
714 |
|
|
ENDIF |
715 |
|
|
if (fl_cor_ebil .LE. 0) then |
716 |
|
|
! the following line should not be here. Indeed, if zifl is modified |
717 |
|
|
! now, zifl(i)*zmelt is no more the amount of ice that has melt |
718 |
|
|
! and therefore the change in temperature computed below is wrong |
719 |
|
|
zifl(i)=zifl(i)*(1.-zmelt) |
720 |
|
|
end if |
721 |
|
|
! Chaleur latente de fusion |
722 |
|
|
if (fl_cor_ebil .GT. 0) then ! avec correction thermalisation des precips |
723 |
|
|
zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) & |
724 |
|
|
*RLMLT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) |
725 |
|
|
else ! sans correction thermalisation des precips |
726 |
|
|
zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) & |
727 |
|
|
*RLMLT/RCPD/(1.0+RVTMP2*zq(i)) |
728 |
|
|
end if |
729 |
|
|
if (fl_cor_ebil .GT. 0) then ! correction bug, deplacement ligne precedente |
730 |
|
|
!<LTP |
731 |
|
|
IF (iflag_evap_prec==4) THEN |
732 |
|
|
ziflclr(i)=ziflclr(i)*(1.-zmelt) |
733 |
|
|
ziflcld(i)=ziflcld(i)*(1.-zmelt) |
734 |
|
|
zifl(i)=ziflclr(i)+ziflcld(i) |
735 |
|
|
!>LTP |
736 |
|
|
ELSE |
737 |
|
|
zifl(i)=zifl(i)*(1.-zmelt) |
738 |
|
|
ENDIF |
739 |
|
|
end if |
740 |
|
|
|
741 |
|
|
ELSE |
742 |
|
|
! Si on n'a plus de pluies, on reinitialise a 0 la farcion |
743 |
|
|
! sous nuageuse utilisee pour la pluie. |
744 |
|
|
znebprecip(i)=0. |
745 |
|
|
ENDIF ! (zrfl(i)+zifl(i).GT.0.) |
746 |
|
|
ENDDO |
747 |
|
|
|
748 |
|
|
ENDIF ! (.NOT. ice_thermo) |
749 |
|
|
|
750 |
|
|
! ---------------------------------------------------------------- |
751 |
|
|
! Fin evaporation de la precipitation |
752 |
|
|
! ---------------------------------------------------------------- |
753 |
|
|
ENDIF ! (iflag_evap_prec>=1) |
754 |
|
|
! |
755 |
|
|
! Calculer Qs et L/Cp*dQs/dT: |
756 |
|
|
! |
757 |
|
|
IF (thermcep) THEN |
758 |
|
|
DO i = 1, klon |
759 |
|
|
zdelta = MAX(0.,SIGN(1.,RTT-zt(i))) |
760 |
|
|
zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta |
761 |
|
|
if (fl_cor_ebil .GT. 0) then ! nouveau |
762 |
|
|
zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) |
763 |
|
|
else |
764 |
|
|
zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i)) |
765 |
|
|
endif |
766 |
|
|
zqs(i) = R2ES*FOEEW(zt(i),zdelta)/pplay(i,k) |
767 |
|
|
zqs(i) = MIN(0.5,zqs(i)) |
768 |
|
|
zcor = 1./(1.-RETV*zqs(i)) |
769 |
|
|
zqs(i) = zqs(i)*zcor |
770 |
|
|
zdqs(i) = FOEDE(zt(i),zdelta,zcvm5,zqs(i),zcor) |
771 |
|
|
zdqsdT_raw(i) = zdqs(i)* & |
772 |
|
|
& RCPD*(1.0+RVTMP2*zq(i)) / (RLVTT*(1.-zdelta) + RLSTT*zdelta) |
773 |
|
|
ENDDO |
774 |
|
|
ELSE |
775 |
|
|
DO i = 1, klon |
776 |
|
|
IF (zt(i).LT.t_coup) THEN |
777 |
|
|
zqs(i) = qsats(zt(i))/pplay(i,k) |
778 |
|
|
zdqs(i) = dqsats(zt(i),zqs(i)) |
779 |
|
|
ELSE |
780 |
|
|
zqs(i) = qsatl(zt(i))/pplay(i,k) |
781 |
|
|
zdqs(i) = dqsatl(zt(i),zqs(i)) |
782 |
|
|
ENDIF |
783 |
|
|
ENDDO |
784 |
|
|
ENDIF |
785 |
|
|
! |
786 |
|
|
! Determiner la condensation partielle et calculer la quantite |
787 |
|
|
! de l'eau condensee: |
788 |
|
|
! |
789 |
|
|
!verification de la valeur de iflag_fisrtilp_qsat pour iflag_ice_thermo=1 |
790 |
|
|
! if ((iflag_ice_thermo.eq.1).and.(iflag_fisrtilp_qsat.ne.0)) then |
791 |
|
|
! write(*,*) " iflag_ice_thermo==1 requires iflag_fisrtilp_qsat==0", & |
792 |
|
|
! " but iflag_fisrtilp_qsat=",iflag_fisrtilp_qsat, ". Might as well stop here." |
793 |
|
|
! stop |
794 |
|
|
! endif |
795 |
|
|
|
796 |
|
|
! ---------------------------------------------------------------- |
797 |
|
|
! P2> Formation du nuage |
798 |
|
|
! ---------------------------------------------------------------- |
799 |
|
|
! Variables calculees: |
800 |
|
|
! rneb : fraction nuageuse |
801 |
|
|
! zcond : eau condensee moyenne dans la maille. |
802 |
|
|
! rhcl: humidite relative ciel-clair |
803 |
|
|
! zt : temperature de la maille |
804 |
|
|
! ---------------------------------------------------------------- |
805 |
|
|
! |
806 |
|
|
IF (cpartiel) THEN |
807 |
|
|
! ------------------------- |
808 |
|
|
! P2.A> Nuage fractionnaire |
809 |
|
|
! ------------------------- |
810 |
|
|
! |
811 |
|
|
! Calcul de l'eau condensee et de la fraction nuageuse et de l'eau |
812 |
|
|
! nuageuse a partir des PDF de Sandrine Bony. |
813 |
|
|
! rneb : fraction nuageuse |
814 |
|
|
! zqn : eau totale dans le nuage |
815 |
|
|
! zcond : eau condensee moyenne dans la maille. |
816 |
|
|
! on prend en compte le réchauffement qui diminue la partie |
817 |
|
|
! condensee |
818 |
|
|
! |
819 |
|
|
! Version avec les raqts |
820 |
|
|
|
821 |
|
|
! ---------------------------------------------------------------- |
822 |
|
|
! P2.A.0> Calcul des grandeurs nuageuses une pdf en creneau |
823 |
|
|
! ---------------------------------------------------------------- |
824 |
|
|
if (iflag_pdf.eq.0) then |
825 |
|
|
|
826 |
|
|
! version creneau de (Li, 1998) |
827 |
|
|
do i=1,klon |
828 |
|
|
zdelq = min(ratqs(i,k),0.99) * zq(i) |
829 |
|
|
rneb(i,k) = (zq(i)+zdelq-zqs(i)) / (2.0*zdelq) |
830 |
|
|
zqn(i) = (zq(i)+zdelq+zqs(i))/2.0 |
831 |
|
|
enddo |
832 |
|
|
|
833 |
|
|
else ! if (iflag_pdf.eq.0) |
834 |
|
|
! ---------------------------------------------------------------- |
835 |
|
|
! P2.A.1> Avec les nouvelles PDFs, calcul des grandeurs nuageuses pour |
836 |
|
|
! les valeurs de T et Q initiales |
837 |
|
|
! ---------------------------------------------------------------- |
838 |
|
|
do i=1,klon |
839 |
|
|
if(zq(i).lt.1.e-15) then |
840 |
|
|
ncoreczq=ncoreczq+1 |
841 |
|
|
zq(i)=1.e-15 |
842 |
|
|
endif |
843 |
|
|
enddo |
844 |
|
|
|
845 |
|
|
if (iflag_cld_th>=5) then |
846 |
|
|
|
847 |
|
|
if (iflag_cloudth_vert<=2) then |
848 |
|
|
call cloudth(klon,klev,k,ztv, & |
849 |
|
|
zq,zqta,fraca, & |
850 |
|
|
qcloud,ctot,zpspsk,paprs,pplay,ztla,zthl, & |
851 |
|
|
ratqs,zqs,t) |
852 |
|
|
elseif (iflag_cloudth_vert>=3 .and. iflag_cloudth_vert<=5) then |
853 |
|
|
call cloudth_v3(klon,klev,k,ztv, & |
854 |
|
|
zq,zqta,fraca, & |
855 |
|
|
qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & |
856 |
|
|
ratqs,zqs,t) |
857 |
|
|
!---------------------------------- |
858 |
|
|
!Version these Jean Jouhaud, Decembre 2018 |
859 |
|
|
!---------------------------------- |
860 |
|
|
elseif (iflag_cloudth_vert==6) then |
861 |
|
|
call cloudth_v6(klon,klev,k,ztv, & |
862 |
|
|
zq,zqta,fraca, & |
863 |
|
|
qcloud,ctot,ctot_vol,zpspsk,paprs,pplay,ztla,zthl, & |
864 |
|
|
ratqs,zqs,t) |
865 |
|
|
|
866 |
|
|
endif |
867 |
|
|
do i=1,klon |
868 |
|
|
rneb(i,k)=ctot(i,k) |
869 |
|
|
rneblsvol(i,k)=ctot_vol(i,k) |
870 |
|
|
zqn(i)=qcloud(i) |
871 |
|
|
enddo |
872 |
|
|
|
873 |
|
|
endif |
874 |
|
|
|
875 |
|
|
if (iflag_cld_th <= 4) then |
876 |
|
|
lognormale = .true. |
877 |
|
|
elseif (iflag_cld_th >= 6) then |
878 |
|
|
! lognormale en l'absence des thermiques |
879 |
|
|
lognormale = fraca(:,k) < 1e-10 |
880 |
|
|
else |
881 |
|
|
! Dans le cas iflag_cld_th=5, on prend systématiquement la |
882 |
|
|
! bi-gaussienne |
883 |
|
|
lognormale = .false. |
884 |
|
|
end if |
885 |
|
|
|
886 |
|
|
!CR: variation de qsat avec T en presence de glace ou non |
887 |
|
|
!initialisations |
888 |
|
|
do i=1,klon |
889 |
|
|
DT(i) = 0. |
890 |
|
|
n_i(i)=0 |
891 |
|
|
Tbef(i)=zt(i) |
892 |
|
|
qlbef(i)=0. |
893 |
|
|
enddo |
894 |
|
|
|
895 |
|
|
! ---------------------------------------------------------------- |
896 |
|
|
! P2.A.2> Prise en compte du couplage entre eau condensee et T. |
897 |
|
|
! Calcul des grandeurs nuageuses en tenant compte de l'effet de |
898 |
|
|
! la condensation sur T, et donc sur qsat et sur les grandeurs nuageuses |
899 |
|
|
! qui en dependent. Ce changement de temperature est provisoire, et |
900 |
|
|
! la valeur definitive sera calcule plus tard. |
901 |
|
|
! Variables calculees: |
902 |
|
|
! rneb : nebulosite |
903 |
|
|
! zcond: eau condensee en moyenne dans la maille |
904 |
|
|
! note JLD: si on n'a pas de pdf lognormale, ce qui se passe ne me semble |
905 |
|
|
! pas clair, il n'y a probablement pas de prise en compte de l'effet de |
906 |
|
|
! T sur qsat |
907 |
|
|
! ---------------------------------------------------------------- |
908 |
|
|
|
909 |
|
|
!Boucle iterative: ATTENTION, l'option -1 n'est plus activable ici |
910 |
|
|
if (iflag_fisrtilp_qsat.ge.0) then |
911 |
|
|
! Iteration pour condensation avec variation de qsat(T) |
912 |
|
|
! ----------------------------------------------------- |
913 |
|
|
do iter=1,iflag_fisrtilp_qsat+1 |
914 |
|
|
|
915 |
|
|
do i=1,klon |
916 |
|
|
! do while ((abs(DT(i)).gt.DDT0).or.(n_i(i).eq.0)) |
917 |
|
|
! !! convergence = .true. tant que l'on n'a pas converge !! |
918 |
|
|
! ------------------------------ |
919 |
|
|
convergence(i)=abs(DT(i)).gt.DDT0 |
920 |
|
|
if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then |
921 |
|
|
! si on n'a pas converge |
922 |
|
|
! |
923 |
|
|
! P2.A.2.1> Calcul de la fraction nuageuse et de la quantite d'eau condensee |
924 |
|
|
! --------------------------------------------------------------- |
925 |
|
|
! Variables calculees: |
926 |
|
|
! rneb : nebulosite |
927 |
|
|
! zqn : eau condensee, dans le nuage (in cloud water content) |
928 |
|
|
! zcond: eau condensee en moyenne dans la maille |
929 |
|
|
! rhcl: humidite relative ciel-clair |
930 |
|
|
! |
931 |
|
|
Tbef(i)=Tbef(i)+DT(i) ! nouvelle temperature |
932 |
|
|
if (.not.ice_thermo) then |
933 |
|
|
zdelta = MAX(0.,SIGN(1.,RTT-Tbef(i))) |
934 |
|
|
else |
935 |
|
|
if (iflag_t_glace.eq.0) then |
936 |
|
|
zdelta = MAX(0.,SIGN(1.,t_glace_min_old-Tbef(i))) |
937 |
|
|
else if (iflag_t_glace.ge.1) then |
938 |
|
|
if (iflag_oldbug_fisrtilp.EQ.0) then |
939 |
|
|
zdelta = MAX(0.,SIGN(1.,t_glace_max-Tbef(i))) |
940 |
|
|
else |
941 |
|
|
!avec bug : zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i))) |
942 |
|
|
zdelta = MAX(0.,SIGN(1.,t_glace_min-Tbef(i))) |
943 |
|
|
endif |
944 |
|
|
endif |
945 |
|
|
endif |
946 |
|
|
! Calcul de rneb, qzn et zcond pour les PDF lognormales |
947 |
|
|
zcvm5 = R5LES*RLVTT*(1.-zdelta) + R5IES*RLSTT*zdelta |
948 |
|
|
if (fl_cor_ebil .GT. 0) then |
949 |
|
|
zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) |
950 |
|
|
else |
951 |
|
|
zcvm5 = zcvm5 /RCPD/(1.0+RVTMP2*zq(i)) |
952 |
|
|
end if |
953 |
|
|
zqs(i) = R2ES*FOEEW(Tbef(i),zdelta)/pplay(i,k) |
954 |
|
|
zqs(i) = MIN(0.5,zqs(i)) |
955 |
|
|
zcor = 1./(1.-RETV*zqs(i)) |
956 |
|
|
zqs(i) = zqs(i)*zcor |
957 |
|
|
zdqs(i) = FOEDE(Tbef(i),zdelta,zcvm5,zqs(i),zcor) |
958 |
|
|
zpdf_sig(i)=ratqs(i,k)*zq(i) |
959 |
|
|
zpdf_k(i)=-sqrt(log(1.+(zpdf_sig(i)/zq(i))**2)) |
960 |
|
|
zpdf_delta(i)=log(zq(i)/zqs(i)) |
961 |
|
|
zpdf_a(i)=zpdf_delta(i)/(zpdf_k(i)*sqrt(2.)) |
962 |
|
|
zpdf_b(i)=zpdf_k(i)/(2.*sqrt(2.)) |
963 |
|
|
zpdf_e1(i)=zpdf_a(i)-zpdf_b(i) |
964 |
|
|
zpdf_e1(i)=sign(min(abs(zpdf_e1(i)),5.),zpdf_e1(i)) |
965 |
|
|
zpdf_e1(i)=1.-erf(zpdf_e1(i)) |
966 |
|
|
zpdf_e2(i)=zpdf_a(i)+zpdf_b(i) |
967 |
|
|
zpdf_e2(i)=sign(min(abs(zpdf_e2(i)),5.),zpdf_e2(i)) |
968 |
|
|
zpdf_e2(i)=1.-erf(zpdf_e2(i)) |
969 |
|
|
|
970 |
|
|
if (zpdf_e1(i).lt.1.e-10) then |
971 |
|
|
rneb(i,k)=0. |
972 |
|
|
zqn(i)=zqs(i) |
973 |
|
|
else |
974 |
|
|
rneb(i,k)=0.5*zpdf_e1(i) |
975 |
|
|
zqn(i)=zq(i)*zpdf_e2(i)/zpdf_e1(i) |
976 |
|
|
endif |
977 |
|
|
|
978 |
|
|
! If vertical heterogeneity, change fraction by volume as well |
979 |
|
|
if (iflag_cloudth_vert>=3) then |
980 |
|
|
ctot_vol(i,k)=rneb(i,k) |
981 |
|
|
rneblsvol(i,k)=ctot_vol(i,k) |
982 |
|
|
endif |
983 |
|
|
|
984 |
|
|
endif !convergence |
985 |
|
|
|
986 |
|
|
enddo ! boucle en i |
987 |
|
|
|
988 |
|
|
! P2.A.2.2> Calcul APPROCHE de la variation de temperature DT |
989 |
|
|
! due a la condensation. |
990 |
|
|
! --------------------------------------------------------------- |
991 |
|
|
! Variables calculees: |
992 |
|
|
! DT : variation de temperature due a la condensation |
993 |
|
|
|
994 |
|
|
if (.not. ice_thermo) then |
995 |
|
|
! -------------------------- |
996 |
|
|
do i=1,klon |
997 |
|
|
if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then |
998 |
|
|
|
999 |
|
|
qlbef(i)=max(0.,zqn(i)-zqs(i)) |
1000 |
|
|
if (fl_cor_ebil .GT. 0) then |
1001 |
|
|
num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i) |
1002 |
|
|
else |
1003 |
|
|
num=-Tbef(i)+zt(i)+rneb(i,k)*RLVTT/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i) |
1004 |
|
|
end if |
1005 |
|
|
denom=1.+rneb(i,k)*zdqs(i) |
1006 |
|
|
DT(i)=num/denom |
1007 |
|
|
n_i(i)=n_i(i)+1 |
1008 |
|
|
endif |
1009 |
|
|
enddo |
1010 |
|
|
|
1011 |
|
|
else ! if (.not. ice_thermo) |
1012 |
|
|
! -------------------------- |
1013 |
|
|
if (iflag_t_glace.ge.1) then |
1014 |
|
|
CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:)) |
1015 |
|
|
endif |
1016 |
|
|
|
1017 |
|
|
do i=1,klon |
1018 |
|
|
if ((convergence(i).or.(n_i(i).eq.0)).and.lognormale(i)) then |
1019 |
|
|
|
1020 |
|
|
if (iflag_t_glace.eq.0) then |
1021 |
|
|
zfice(i) = 1.0 - (Tbef(i)-t_glace_min_old) / (RTT-t_glace_min_old) |
1022 |
|
|
zfice(i) = MIN(MAX(zfice(i),0.0),1.0) |
1023 |
|
|
zfice(i) = zfice(i)**exposant_glace_old |
1024 |
|
|
dzfice(i)= exposant_glace_old * zfice(i)**(exposant_glace_old-1) & |
1025 |
|
|
& / (t_glace_min_old - RTT) |
1026 |
|
|
endif |
1027 |
|
|
|
1028 |
|
|
if (iflag_t_glace.ge.1.and.zfice(i)>0.) then |
1029 |
|
|
dzfice(i)= exposant_glace * zfice(i)**(exposant_glace-1) & |
1030 |
|
|
& / (t_glace_min - t_glace_max) |
1031 |
|
|
endif |
1032 |
|
|
|
1033 |
|
|
if ((zfice(i).eq.0).or.(zfice(i).eq.1)) then |
1034 |
|
|
dzfice(i)=0. |
1035 |
|
|
endif |
1036 |
|
|
|
1037 |
|
|
if (zfice(i).lt.1) then |
1038 |
|
|
cste=RLVTT |
1039 |
|
|
else |
1040 |
|
|
cste=RLSTT |
1041 |
|
|
endif |
1042 |
|
|
|
1043 |
|
|
qlbef(i)=max(0.,zqn(i)-zqs(i)) |
1044 |
|
|
if (fl_cor_ebil .GT. 0) then |
1045 |
|
|
num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT & |
1046 |
|
|
& +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*qlbef(i) |
1047 |
|
|
denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) & |
1048 |
|
|
-(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)))*rneb(i,k) & |
1049 |
|
|
& *qlbef(i)*dzfice(i) |
1050 |
|
|
else |
1051 |
|
|
num = -Tbef(i)+zt(i)+rneb(i,k)*((1-zfice(i))*RLVTT & |
1052 |
|
|
& +zfice(i)*RLSTT)/RCPD/(1.0+RVTMP2*zq(i))*qlbef(i) |
1053 |
|
|
denom = 1.+rneb(i,k)*((1-zfice(i))*RLVTT+zfice(i)*RLSTT)/cste*zdqs(i) & |
1054 |
|
|
-(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i))*rneb(i,k)*qlbef(i)*dzfice(i) |
1055 |
|
|
end if |
1056 |
|
|
DT(i)=num/denom |
1057 |
|
|
n_i(i)=n_i(i)+1 |
1058 |
|
|
|
1059 |
|
|
endif ! fin convergence |
1060 |
|
|
enddo ! fin boucle i |
1061 |
|
|
|
1062 |
|
|
endif !ice_thermo |
1063 |
|
|
|
1064 |
|
|
enddo ! iter=1,iflag_fisrtilp_qsat+1 |
1065 |
|
|
! Fin d'iteration pour condensation avec variation de qsat(T) |
1066 |
|
|
! ----------------------------------------------------------- |
1067 |
|
|
endif ! if (iflag_fisrtilp_qsat.ge.0) |
1068 |
|
|
! ---------------------------------------------------------------- |
1069 |
|
|
! Fin de P2.A.2> la prise en compte du couplage entre eau condensee et T |
1070 |
|
|
! ---------------------------------------------------------------- |
1071 |
|
|
|
1072 |
|
|
endif ! iflag_pdf |
1073 |
|
|
|
1074 |
|
|
! if (iflag_fisrtilp_qsat.eq.-1) then |
1075 |
|
|
!------------------------------------------ |
1076 |
|
|
!CR: ATTENTION option fausse mais a existe: |
1077 |
|
|
! pour la re-activer, prendre iflag_fisrtilp_qsat=0 et |
1078 |
|
|
! activer les lignes suivantes: |
1079 |
|
|
IF (1.eq.0) THEN |
1080 |
|
|
DO i=1,klon |
1081 |
|
|
IF (rneb(i,k) .LE. 0.0) THEN |
1082 |
|
|
zqn(i) = 0.0 |
1083 |
|
|
rneb(i,k) = 0.0 |
1084 |
|
|
zcond(i) = 0.0 |
1085 |
|
|
rhcl(i,k)=zq(i)/zqs(i) |
1086 |
|
|
ELSE IF (rneb(i,k) .GE. 1.0) THEN |
1087 |
|
|
zqn(i) = zq(i) |
1088 |
|
|
rneb(i,k) = 1.0 |
1089 |
|
|
zcond(i) = MAX(0.0,zqn(i)-zqs(i))/(1+zdqs(i)) |
1090 |
|
|
rhcl(i,k)=1.0 |
1091 |
|
|
ELSE |
1092 |
|
|
zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k)/(1+zdqs(i)) |
1093 |
|
|
rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i) |
1094 |
|
|
ENDIF |
1095 |
|
|
ENDDO |
1096 |
|
|
ENDIF |
1097 |
|
|
!------------------------------------------ |
1098 |
|
|
|
1099 |
|
|
! ELSE |
1100 |
|
|
! ---------------------------------------------------------------- |
1101 |
|
|
! P2.A.3> Calcul des valeures finales associees a la formation des nuages |
1102 |
|
|
! Variables calculees: |
1103 |
|
|
! rneb : nebulosite |
1104 |
|
|
! zcond: eau condensee en moyenne dans la maille |
1105 |
|
|
! zq : eau vapeur dans la maille |
1106 |
|
|
! zt : temperature de la maille |
1107 |
|
|
! rhcl: humidite relative ciel-clair |
1108 |
|
|
! ---------------------------------------------------------------- |
1109 |
|
|
! |
1110 |
|
|
! Bornage de l'eau in-cloud (zqn) et de la fraction nuageuse (rneb) |
1111 |
|
|
! Calcule de l'eau condensee moyenne dans la maille (zcond), |
1112 |
|
|
! et de l'humidite relative ciel-clair (rhcl) |
1113 |
|
|
DO i=1,klon |
1114 |
|
|
IF (rneb(i,k) .LE. 0.0) THEN |
1115 |
|
|
zqn(i) = 0.0 |
1116 |
|
|
rneb(i,k) = 0.0 |
1117 |
|
|
zcond(i) = 0.0 |
1118 |
|
|
rhcl(i,k)=zq(i)/zqs(i) |
1119 |
|
|
ELSE IF (rneb(i,k) .GE. 1.0) THEN |
1120 |
|
|
zqn(i) = zq(i) |
1121 |
|
|
rneb(i,k) = 1.0 |
1122 |
|
|
zcond(i) = MAX(0.0,zqn(i)-zqs(i)) |
1123 |
|
|
rhcl(i,k)=1.0 |
1124 |
|
|
ELSE |
1125 |
|
|
zcond(i) = MAX(0.0,zqn(i)-zqs(i))*rneb(i,k) |
1126 |
|
|
rhcl(i,k)=(zqs(i)+zq(i)-zdelq)/2./zqs(i) |
1127 |
|
|
ENDIF |
1128 |
|
|
ENDDO |
1129 |
|
|
|
1130 |
|
|
|
1131 |
|
|
! If vertical heterogeneity, change fraction by volume as well |
1132 |
|
|
if (iflag_cloudth_vert>=3) then |
1133 |
|
|
ctot_vol(1:klon,k)=min(max(ctot_vol(1:klon,k),0.),1.) |
1134 |
|
|
rneblsvol(1:klon,k)=ctot_vol(1:klon,k) |
1135 |
|
|
endif |
1136 |
|
|
|
1137 |
|
|
! ENDIF |
1138 |
|
|
|
1139 |
|
|
ELSE ! de IF (cpartiel) |
1140 |
|
|
! ------------------------- |
1141 |
|
|
! P2.B> Nuage "tout ou rien" |
1142 |
|
|
! ------------------------- |
1143 |
|
|
! note JLD: attention, rhcl non calcule. Ca peut avoir des consequences? |
1144 |
|
|
DO i = 1, klon |
1145 |
|
|
IF (zq(i).GT.zqs(i)) THEN |
1146 |
|
|
rneb(i,k) = 1.0 |
1147 |
|
|
ELSE |
1148 |
|
|
rneb(i,k) = 0.0 |
1149 |
|
|
ENDIF |
1150 |
|
|
zcond(i) = MAX(0.0,zq(i)-zqs(i))/(1.+zdqs(i)) |
1151 |
|
|
ENDDO |
1152 |
|
|
ENDIF ! de IF (cpartiel) |
1153 |
|
|
! |
1154 |
|
|
! Mise a jour vapeur d'eau |
1155 |
|
|
! ------------------------- |
1156 |
|
|
DO i = 1, klon |
1157 |
|
|
zq(i) = zq(i) - zcond(i) |
1158 |
|
|
! zt(i) = zt(i) + zcond(i) * RLVTT/RCPD |
1159 |
|
|
ENDDO |
1160 |
|
|
!AJ< |
1161 |
|
|
! ------------------------------------ |
1162 |
|
|
! P2.C> Prise en compte de la Chaleur latente apres formation nuage |
1163 |
|
|
! ------------------------------------- |
1164 |
|
|
! Variable calcule: |
1165 |
|
|
! zt : temperature de la maille |
1166 |
|
|
! |
1167 |
|
|
IF (.NOT. ice_thermo) THEN |
1168 |
|
|
if (iflag_fisrtilp_qsat.lt.1) then |
1169 |
|
|
DO i = 1, klon |
1170 |
|
|
zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) |
1171 |
|
|
ENDDO |
1172 |
|
|
else if (iflag_fisrtilp_qsat.gt.0) then |
1173 |
|
|
DO i= 1, klon |
1174 |
|
|
if (fl_cor_ebil .GT. 0) then |
1175 |
|
|
zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) |
1176 |
|
|
else |
1177 |
|
|
zt(i) = zt(i) + zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) |
1178 |
|
|
end if |
1179 |
|
|
ENDDO |
1180 |
|
|
endif |
1181 |
|
|
ELSE |
1182 |
|
|
if (iflag_t_glace.ge.1) then |
1183 |
|
|
CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:)) |
1184 |
|
|
endif |
1185 |
|
|
if (iflag_fisrtilp_qsat.lt.1) then |
1186 |
|
|
DO i = 1, klon |
1187 |
|
|
! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod |
1188 |
|
|
! zfice(i) = icefrac_lsc(zt(i), t_glace_min, & |
1189 |
|
|
! t_glace_max, exposant_glace) |
1190 |
|
|
if (iflag_t_glace.eq.0) then |
1191 |
|
|
zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old) |
1192 |
|
|
zfice(i) = MIN(MAX(zfice(i),0.0),1.0) |
1193 |
|
|
zfice(i) = zfice(i)**exposant_glace_old |
1194 |
|
|
endif |
1195 |
|
|
zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*zq(i)) & |
1196 |
|
|
+zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*zq(i)) |
1197 |
|
|
ENDDO |
1198 |
|
|
else |
1199 |
|
|
DO i=1, klon |
1200 |
|
|
! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod |
1201 |
|
|
! zfice(i) = icefrac_lsc(zt(i), t_glace_min, & |
1202 |
|
|
! t_glace_max, exposant_glace) |
1203 |
|
|
if (iflag_t_glace.eq.0) then |
1204 |
|
|
zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (RTT-t_glace_min_old) |
1205 |
|
|
zfice(i) = MIN(MAX(zfice(i),0.0),1.0) |
1206 |
|
|
zfice(i) = zfice(i)**exposant_glace_old |
1207 |
|
|
endif |
1208 |
|
|
if (fl_cor_ebil .GT. 0) then |
1209 |
|
|
zt(i) = zt(i) + (1.-zfice(i))*zcond(i) & |
1210 |
|
|
& * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) & |
1211 |
|
|
+zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) |
1212 |
|
|
else |
1213 |
|
|
zt(i) = zt(i) + (1.-zfice(i))*zcond(i) * RLVTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) & |
1214 |
|
|
+zfice(i)*zcond(i) * RLSTT/RCPD/(1.0+RVTMP2*(zq(i)+zcond(i))) |
1215 |
|
|
end if |
1216 |
|
|
ENDDO |
1217 |
|
|
endif |
1218 |
|
|
! print*,zt(i),zrfl(i),zifl(i),'temp1' |
1219 |
|
|
ENDIF |
1220 |
|
|
!>AJ |
1221 |
|
|
|
1222 |
|
|
! ---------------------------------------------------------------- |
1223 |
|
|
! P3> Formation des precipitations |
1224 |
|
|
! ---------------------------------------------------------------- |
1225 |
|
|
! |
1226 |
|
|
! Partager l'eau condensee en precipitation et eau liquide nuageuse |
1227 |
|
|
! |
1228 |
|
|
|
1229 |
|
|
!<LTP |
1230 |
|
|
|
1231 |
|
|
IF (iflag_evap_prec==4) THEN |
1232 |
|
|
!Partitionnement des precipitations venant du dessus en précipitations nuageuses |
1233 |
|
|
!et précipitations ciel clair |
1234 |
|
|
|
1235 |
|
|
!0) Calculate tot_zneb, la fraction nuageuse totale au-dessus du nuage |
1236 |
|
|
!en supposant un recouvrement maximum aléatoire (voir Jakob and Klein, 2000) |
1237 |
|
|
|
1238 |
|
|
DO i=1, klon |
1239 |
|
|
tot_znebn(i) = 1 - (1-tot_zneb(i))*(1 - max(rneb(i,k),zneb(i))) & |
1240 |
|
|
/(1-min(zneb(i),1-smallestreal)) |
1241 |
|
|
d_tot_zneb(i) = tot_znebn(i) - tot_zneb(i) |
1242 |
|
|
tot_zneb(i) = tot_znebn(i) |
1243 |
|
|
|
1244 |
|
|
|
1245 |
|
|
!1) Cloudy to clear air |
1246 |
|
|
d_znebprecip_cld_clr(i) = znebprecipcld(i) - min(rneb(i,k),znebprecipcld(i)) |
1247 |
|
|
IF (znebprecipcld(i) .GT. 0) THEN |
1248 |
|
|
d_zrfl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*zrflcld(i) |
1249 |
|
|
d_zifl_cld_clr(i) = d_znebprecip_cld_clr(i)/znebprecipcld(i)*ziflcld(i) |
1250 |
|
|
ELSE |
1251 |
|
|
d_zrfl_cld_clr(i) = 0. |
1252 |
|
|
d_zifl_cld_clr(i) = 0. |
1253 |
|
|
ENDIF |
1254 |
|
|
|
1255 |
|
|
!2) Clear to cloudy air |
1256 |
|
|
d_znebprecip_clr_cld(i) = max(0., min(znebprecipclr(i), rneb(i,k) & |
1257 |
|
|
- d_tot_zneb(i) - zneb(i))) |
1258 |
|
|
IF (znebprecipclr(i) .GT. 0) THEN |
1259 |
|
|
d_zrfl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*zrflclr(i) |
1260 |
|
|
d_zifl_clr_cld(i) = d_znebprecip_clr_cld(i)/znebprecipclr(i)*ziflclr(i) |
1261 |
|
|
ELSE |
1262 |
|
|
d_zrfl_clr_cld(i) = 0. |
1263 |
|
|
d_zifl_clr_cld(i) = 0. |
1264 |
|
|
ENDIF |
1265 |
|
|
|
1266 |
|
|
!Update variables |
1267 |
|
|
znebprecipcld(i) = znebprecipcld(i) + d_znebprecip_clr_cld(i) - d_znebprecip_cld_clr(i) |
1268 |
|
|
znebprecipclr(i) = znebprecipclr(i) + d_znebprecip_cld_clr(i) - d_znebprecip_clr_cld(i) |
1269 |
|
|
zrflcld(i) = zrflcld(i) + d_zrfl_clr_cld(i) - d_zrfl_cld_clr(i) |
1270 |
|
|
ziflcld(i) = ziflcld(i) + d_zifl_clr_cld(i) - d_zifl_cld_clr(i) |
1271 |
|
|
zrflclr(i) = zrflclr(i) + d_zrfl_cld_clr(i) - d_zrfl_clr_cld(i) |
1272 |
|
|
ziflclr(i) = ziflclr(i) + d_zifl_cld_clr(i) - d_zifl_clr_cld(i) |
1273 |
|
|
|
1274 |
|
|
ENDDO |
1275 |
|
|
ENDIF |
1276 |
|
|
|
1277 |
|
|
!>LTP |
1278 |
|
|
|
1279 |
|
|
|
1280 |
|
|
|
1281 |
|
|
! Initialisation de zoliq (eau condensee moyenne dans la maille) |
1282 |
|
|
DO i = 1, klon |
1283 |
|
|
IF (rneb(i,k).GT.0.0) THEN |
1284 |
|
|
zoliq(i) = zcond(i) |
1285 |
|
|
zrho(i) = pplay(i,k) / zt(i) / RD |
1286 |
|
|
zdz(i) = (paprs(i,k)-paprs(i,k+1)) / (zrho(i)*RG) |
1287 |
|
|
ENDIF |
1288 |
|
|
ENDDO |
1289 |
|
|
!AJ< |
1290 |
|
|
IF (.NOT. ice_thermo) THEN |
1291 |
|
|
IF (iflag_t_glace.EQ.0) THEN |
1292 |
|
|
DO i = 1, klon |
1293 |
|
|
IF (rneb(i,k).GT.0.0) THEN |
1294 |
|
|
zfice(i) = 1.0 - (zt(i)-t_glace_min_old) / (273.13-t_glace_min_old) |
1295 |
|
|
zfice(i) = MIN(MAX(zfice(i),0.0),1.0) |
1296 |
|
|
zfice(i) = zfice(i)**exposant_glace_old |
1297 |
|
|
! zfice(i) = zfice(i)**nexpo |
1298 |
|
|
!! zfice(i)=0. |
1299 |
|
|
ENDIF |
1300 |
|
|
ENDDO |
1301 |
|
|
ELSE ! of IF (iflag_t_glace.EQ.0) |
1302 |
|
|
CALL icefrac_lsc(klon,zt(:),pplay(:,k)/paprs(:,1),zfice(:)) |
1303 |
|
|
! DO i = 1, klon |
1304 |
|
|
! IF (rneb(i,k).GT.0.0) THEN |
1305 |
|
|
! JBM: icefrac_lsc is now a function contained in icefrac_lsc_mod |
1306 |
|
|
! zfice(i) = icefrac_lsc(zt(i), t_glace_min, & |
1307 |
|
|
! t_glace_max, exposant_glace) |
1308 |
|
|
! ENDIF |
1309 |
|
|
! ENDDO |
1310 |
|
|
ENDIF |
1311 |
|
|
ENDIF |
1312 |
|
|
|
1313 |
|
|
! Calcul de radliq (eau condensee pour le rayonnement) |
1314 |
|
|
! Iteration pour realiser une moyenne de l'eau nuageuse lors de la precip |
1315 |
|
|
! Remarque: ce n'est donc pas l'eau restante en fin de precip mais une |
1316 |
|
|
! eau moyenne restante dans le nuage sur la duree du pas de temps qui est |
1317 |
|
|
! transmise au rayonnement; |
1318 |
|
|
! ---------------------------------------------------------------- |
1319 |
|
|
DO i = 1, klon |
1320 |
|
|
IF (rneb(i,k).GT.0.0) THEN |
1321 |
|
|
zneb(i) = MAX(rneb(i,k), seuil_neb) |
1322 |
|
|
! zt(i) = zt(i)+zcond(i)*zfice(i)*RLMLT/RCPD/(1.0+RVTMP2*zq(i)) |
1323 |
|
|
! print*,zt(i),'fractionglace' |
1324 |
|
|
!>AJ |
1325 |
|
|
radliq(i,k) = zoliq(i)/REAL(ninter+1) |
1326 |
|
|
ENDIF |
1327 |
|
|
ENDDO |
1328 |
|
|
! |
1329 |
|
|
DO n = 1, ninter |
1330 |
|
|
DO i = 1, klon |
1331 |
|
|
IF (rneb(i,k).GT.0.0) THEN |
1332 |
|
|
zrhol(i) = zrho(i) * zoliq(i) / zneb(i) |
1333 |
|
|
! Initialization of zpluie and zice: |
1334 |
|
|
zpluie=0 |
1335 |
|
|
zice=0 |
1336 |
|
|
IF (zneb(i).EQ.seuil_neb) THEN |
1337 |
|
|
ztot = 0.0 |
1338 |
|
|
ELSE |
1339 |
|
|
! quantite d'eau a eliminer: zchau (Sundqvist, 1978) |
1340 |
|
|
! meme chose pour la glace: zfroi (Zender & Kiehl, 1997) |
1341 |
|
|
if (ptconv(i,k)) then |
1342 |
|
|
zcl =cld_lc_con |
1343 |
|
|
zct =1./cld_tau_con |
1344 |
|
|
zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) & |
1345 |
|
|
*fallvc(zrhol(i)) * zfice(i) |
1346 |
|
|
else |
1347 |
|
|
zcl =cld_lc_lsc |
1348 |
|
|
zct =1./cld_tau_lsc |
1349 |
|
|
zfroi = dtime/REAL(ninter)/zdz(i)*zoliq(i) & |
1350 |
|
|
*fallvs(zrhol(i)) * zfice(i) |
1351 |
|
|
endif |
1352 |
|
|
|
1353 |
|
|
! si l'heterogeneite verticale est active, on utilise |
1354 |
|
|
! la fraction volumique "vraie" plutot que la fraction |
1355 |
|
|
! surfacique modifiee, qui est plus grande et reduit |
1356 |
|
|
! sinon l'eau in-cloud de facon artificielle |
1357 |
|
|
if ((iflag_cloudth_vert>=3).AND.(iflag_rain_incloud_vol==1)) then |
1358 |
|
|
zchau = zct *dtime/REAL(ninter) * zoliq(i) & |
1359 |
|
|
*(1.0-EXP(-(zoliq(i)/ctot_vol(i,k)/zcl )**2)) *(1.-zfice(i)) |
1360 |
|
|
else |
1361 |
|
|
zchau = zct *dtime/REAL(ninter) * zoliq(i) & |
1362 |
|
|
*(1.0-EXP(-(zoliq(i)/zneb(i)/zcl )**2)) *(1.-zfice(i)) |
1363 |
|
|
endif |
1364 |
|
|
!AJ< |
1365 |
|
|
IF (.NOT. ice_thermo) THEN |
1366 |
|
|
ztot = zchau + zfroi |
1367 |
|
|
ELSE |
1368 |
|
|
zpluie = MIN(MAX(zchau,0.0),zoliq(i)*(1.-zfice(i))) |
1369 |
|
|
zice = MIN(MAX(zfroi,0.0),zoliq(i)*zfice(i)) |
1370 |
|
|
ztot = zpluie + zice |
1371 |
|
|
ENDIF |
1372 |
|
|
!>AJ |
1373 |
|
|
ztot = MAX(ztot ,0.0) |
1374 |
|
|
ENDIF |
1375 |
|
|
ztot = MIN(ztot,zoliq(i)) |
1376 |
|
|
!AJ< |
1377 |
|
|
! zoliqp = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie , 0.0) |
1378 |
|
|
! zoliqi = MAX(zoliq(i)*zfice(i)-1.*zice , 0.0) |
1379 |
|
|
!JLD : les 2 variables zoliqp et zoliqi crorresponent a des pseudo precip |
1380 |
|
|
! temporaires et ne doivent pas etre calcule (alors qu'elles le sont |
1381 |
|
|
! si iflag_bergeron <> 2 |
1382 |
|
|
! A SUPPRIMER A TERME |
1383 |
|
|
zoliqp(i) = MAX(zoliq(i)*(1.-zfice(i))-1.*zpluie , 0.0) |
1384 |
|
|
zoliqi(i) = MAX(zoliq(i)*zfice(i)-1.*zice , 0.0) |
1385 |
|
|
zoliq(i) = MAX(zoliq(i)-ztot , 0.0) |
1386 |
|
|
!>AJ |
1387 |
|
|
radliq(i,k) = radliq(i,k) + zoliq(i)/REAL(ninter+1) |
1388 |
|
|
ENDIF |
1389 |
|
|
ENDDO ! i = 1,klon |
1390 |
|
|
ENDDO ! n = 1,ninter |
1391 |
|
|
|
1392 |
|
|
! ---------------------------------------------------------------- |
1393 |
|
|
! |
1394 |
|
|
IF (.NOT. ice_thermo) THEN |
1395 |
|
|
DO i = 1, klon |
1396 |
|
|
IF (rneb(i,k).GT.0.0) THEN |
1397 |
|
|
d_ql(i,k) = zoliq(i) |
1398 |
|
|
zrfl(i) = zrfl(i)+ MAX(zcond(i)-zoliq(i),0.0) & |
1399 |
|
|
* (paprs(i,k)-paprs(i,k+1))/(RG*dtime) |
1400 |
|
|
ENDIF |
1401 |
|
|
ENDDO |
1402 |
|
|
ELSE |
1403 |
|
|
! |
1404 |
|
|
!CR&JYG< |
1405 |
|
|
! On prend en compte l'effet Bergeron dans les flux de precipitation : |
1406 |
|
|
! Si T < 0 C, alors les precipitations liquides sont converties en glace, ce qui |
1407 |
|
|
! provoque un accroissement de temperature DeltaT. L'effet de DeltaT sur le condensat |
1408 |
|
|
! et les precipitations est grossierement pris en compte en linearisant les equations |
1409 |
|
|
! et en approximant le processus de precipitation liquide par un processus a seuil. |
1410 |
|
|
! On fait l'hypothese que le condensat nuageux n'est pas modifié dans cette opération. |
1411 |
|
|
! Le condensat precipitant liquide est supprime (dans la limite DeltaT<273-T). |
1412 |
|
|
! Le condensat precipitant solide est augmente. |
1413 |
|
|
! La vapeur d'eau est augmentee. |
1414 |
|
|
! |
1415 |
|
|
IF ((iflag_bergeron .EQ. 2)) THEN |
1416 |
|
|
DO i = 1, klon |
1417 |
|
|
IF (rneb(i,k) .GT. 0.0) THEN |
1418 |
|
|
zqpreci(i)=(zcond(i)-zoliq(i))*zfice(i) |
1419 |
|
|
zqprecl(i)=(zcond(i)-zoliq(i))*(1.-zfice(i)) |
1420 |
|
|
if (fl_cor_ebil .GT. 0) then |
1421 |
|
|
zcp=RCPD*(1.0+RVTMP2*(zq(i)+zmqc(i)+zcond(i))) |
1422 |
|
|
coef1 = rneb(i,k)*RLSTT/zcp*zdqsdT_raw(i) |
1423 |
|
|
! Calcul de DT si toute les precips liquides congelent |
1424 |
|
|
DeltaT = RLMLT*zqprecl(i) / (zcp*(1.+coef1)) |
1425 |
|
|
! On ne veut pas que T devienne superieur a la temp. de congelation. |
1426 |
|
|
! donc que Delta > RTT-zt(i |
1427 |
|
|
DeltaT = max( min( RTT-zt(i), DeltaT) , 0. ) |
1428 |
|
|
zt(i) = zt(i) + DeltaT |
1429 |
|
|
! Eau vaporisee du fait de l'augmentation de T |
1430 |
|
|
Deltaq = rneb(i,k)*zdqsdT_raw(i)*DeltaT |
1431 |
|
|
! on reajoute cette eau vaporise a la vapeur et on l'enleve des precips |
1432 |
|
|
zq(i) = zq(i) + Deltaq |
1433 |
|
|
! Les 3 max si dessous prtotegent uniquement des erreurs d'arrondies |
1434 |
|
|
zcond(i) = max( zcond(i)- Deltaq, 0. ) |
1435 |
|
|
! precip liquide qui congele ou qui s'evapore |
1436 |
|
|
Deltaqprecl = -zcp/RLMLT*(1.+coef1)*DeltaT |
1437 |
|
|
zqprecl(i) = max( zqprecl(i) + Deltaqprecl, 0. ) |
1438 |
|
|
! bilan eau glacee |
1439 |
|
|
zqpreci(i) = max (zqpreci(i) - Deltaqprecl - Deltaq, 0.) |
1440 |
|
|
else ! if (fl_cor_ebil .GT. 0) |
1441 |
|
|
! ancien calcul |
1442 |
|
|
zcp=RCPD*(1.0+RVTMP2*(zq(i)+zcond(i))) |
1443 |
|
|
coef1 = RLMLT*zdqs(i)/RLVTT |
1444 |
|
|
DeltaT = max( min( RTT-zt(i), RLMLT*zqprecl(i)/zcp/(1.+coef1) ) , 0.) |
1445 |
|
|
zqpreci(i) = zqpreci(i) + zcp/RLMLT*DeltaT |
1446 |
|
|
zqprecl(i) = max( zqprecl(i) - zcp/RLMLT*(1.+coef1)*DeltaT, 0. ) |
1447 |
|
|
zcond(i) = max( zcond(i) - zcp/RLVTT*zdqs(i)*DeltaT, 0. ) |
1448 |
|
|
zq(i) = zq(i) + zcp/RLVTT*zdqs(i)*DeltaT |
1449 |
|
|
zt(i) = zt(i) + DeltaT |
1450 |
|
|
end if ! if (fl_cor_ebil .GT. 0) |
1451 |
|
|
ENDIF ! rneb(i,k) .GT. 0.0 |
1452 |
|
|
ENDDO |
1453 |
|
|
DO i = 1, klon |
1454 |
|
|
IF (rneb(i,k).GT.0.0) THEN |
1455 |
|
|
d_ql(i,k) = (1-zfice(i))*zoliq(i) |
1456 |
|
|
d_qi(i,k) = zfice(i)*zoliq(i) |
1457 |
|
|
!<LTP |
1458 |
|
|
IF (iflag_evap_prec == 4) THEN |
1459 |
|
|
zrflcld(i) = zrflcld(i)+zqprecl(i) & |
1460 |
|
|
*(paprs(i,k)-paprs(i,k+1))/(RG*dtime) |
1461 |
|
|
ziflcld(i) = ziflcld(i)+ zqpreci(i) & |
1462 |
|
|
*(paprs(i,k)-paprs(i,k+1))/(RG*dtime) |
1463 |
|
|
znebprecipcld(i) = rneb(i,k) |
1464 |
|
|
zrfl(i) = zrflcld(i) + zrflclr(i) |
1465 |
|
|
zifl(i) = ziflcld(i) + ziflclr(i) |
1466 |
|
|
!>LTP |
1467 |
|
|
ELSE |
1468 |
|
|
zrfl(i) = zrfl(i)+ zqprecl(i) & |
1469 |
|
|
*(paprs(i,k)-paprs(i,k+1))/(RG*dtime) |
1470 |
|
|
zifl(i) = zifl(i)+ zqpreci(i) & |
1471 |
|
|
*(paprs(i,k)-paprs(i,k+1))/(RG*dtime) |
1472 |
|
|
|
1473 |
|
|
ENDIF !iflag_evap_prec==4 |
1474 |
|
|
|
1475 |
|
|
ENDIF |
1476 |
|
|
ENDDO |
1477 |
|
|
!! |
1478 |
|
|
ELSE ! iflag_bergeron |
1479 |
|
|
!>CR&JYG |
1480 |
|
|
!! |
1481 |
|
|
DO i = 1, klon |
1482 |
|
|
IF (rneb(i,k).GT.0.0) THEN |
1483 |
|
|
!CR on prend en compte la phase glace |
1484 |
|
|
!JLD inutile car on ne passe jamais ici si .not.ice_thermo |
1485 |
|
|
! if (.not.ice_thermo) then |
1486 |
|
|
! d_ql(i,k) = zoliq(i) |
1487 |
|
|
! d_qi(i,k) = 0. |
1488 |
|
|
! else |
1489 |
|
|
d_ql(i,k) = (1-zfice(i))*zoliq(i) |
1490 |
|
|
d_qi(i,k) = zfice(i)*zoliq(i) |
1491 |
|
|
! endif |
1492 |
|
|
!<LTP |
1493 |
|
|
IF (iflag_evap_prec == 4) THEN |
1494 |
|
|
zrflcld(i) = zrflcld(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) & |
1495 |
|
|
*(paprs(i,k)-paprs(i,k+1))/(RG*dtime) |
1496 |
|
|
ziflcld(i) = ziflcld(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) & |
1497 |
|
|
*(paprs(i,k)-paprs(i,k+1))/(RG*dtime) |
1498 |
|
|
znebprecipcld(i) = rneb(i,k) |
1499 |
|
|
zrfl(i) = zrflcld(i) + zrflclr(i) |
1500 |
|
|
zifl(i) = ziflcld(i) + ziflclr(i) |
1501 |
|
|
!>LTP |
1502 |
|
|
ELSE |
1503 |
|
|
!AJ< |
1504 |
|
|
zrfl(i) = zrfl(i)+ MAX(zcond(i)*(1.-zfice(i))-zoliqp(i),0.0) & |
1505 |
|
|
*(paprs(i,k)-paprs(i,k+1))/(RG*dtime) |
1506 |
|
|
zifl(i) = zifl(i)+ MAX(zcond(i)*zfice(i)-zoliqi(i),0.0) & |
1507 |
|
|
*(paprs(i,k)-paprs(i,k+1))/(RG*dtime) |
1508 |
|
|
! zrfl(i) = zrfl(i)+ zpluie & |
1509 |
|
|
! *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) |
1510 |
|
|
! zifl(i) = zifl(i)+ zice & |
1511 |
|
|
! *(paprs(i,k)-paprs(i,k+1))/(RG*dtime) |
1512 |
|
|
ENDIF !iflag_evap_prec == 4 |
1513 |
|
|
|
1514 |
|
|
!CR : on prend en compte l'effet Bergeron dans les flux de precipitation |
1515 |
|
|
IF ((iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15)) THEN |
1516 |
|
|
!<LTP |
1517 |
|
|
IF (iflag_evap_prec == 4) THEN |
1518 |
|
|
zsolid = zrfl(i) |
1519 |
|
|
ziflclr(i) = ziflclr(i) +zrflclr(i) |
1520 |
|
|
ziflcld(i) = ziflcld(i) +zrflcld(i) |
1521 |
|
|
zifl(i) = ziflclr(i)+ziflcld(i) |
1522 |
|
|
zrflcld(i)=0. |
1523 |
|
|
zrflclr(i)=0. |
1524 |
|
|
zrfl(i) = zrflclr(i)+zrflcld(i) |
1525 |
|
|
!>LTP |
1526 |
|
|
ELSE |
1527 |
|
|
zsolid = zrfl(i) |
1528 |
|
|
zifl(i) = zifl(i)+zrfl(i) |
1529 |
|
|
zrfl(i) = 0. |
1530 |
|
|
ENDIF!iflag_evap_prec==4 |
1531 |
|
|
|
1532 |
|
|
if (fl_cor_ebil .GT. 0) then |
1533 |
|
|
zt(i)=zt(i)+zsolid*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) & |
1534 |
|
|
*(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*(zq(i)+zmqc(i))) |
1535 |
|
|
else |
1536 |
|
|
zt(i)=zt(i)+zsolid*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) & |
1537 |
|
|
*(RLSTT-RLVTT)/RCPD/(1.0+RVTMP2*zq(i)) |
1538 |
|
|
end if |
1539 |
|
|
ENDIF ! (iflag_bergeron .EQ. 1) .AND. (zt(i) .LT. 273.15) |
1540 |
|
|
!RC |
1541 |
|
|
|
1542 |
|
|
ENDIF ! rneb(i,k).GT.0.0 |
1543 |
|
|
ENDDO |
1544 |
|
|
|
1545 |
|
|
ENDIF ! iflag_bergeron .EQ. 2 |
1546 |
|
|
ENDIF ! .NOT. ice_thermo |
1547 |
|
|
|
1548 |
|
|
!CR: la fonte est faite au debut |
1549 |
|
|
! IF (ice_thermo) THEN |
1550 |
|
|
! DO i = 1, klon |
1551 |
|
|
! zmelt = ((zt(i)-273.15)/(ztfondue-273.15))**2 |
1552 |
|
|
! zmelt = MIN(MAX(zmelt,0.),1.) |
1553 |
|
|
! zrfl(i)=zrfl(i)+zmelt*zifl(i) |
1554 |
|
|
! zifl(i)=zifl(i)*(1.-zmelt) |
1555 |
|
|
! print*,zt(i),'octavio1' |
1556 |
|
|
! zt(i)=zt(i)-zifl(i)*zmelt*(RG*dtime)/(paprs(i,k)-paprs(i,k+1)) & |
1557 |
|
|
! *RLMLT/RCPD/(1.0+RVTMP2*zq(i)) |
1558 |
|
|
! print*,zt(i),zrfl(i),zifl(i),zmelt,'octavio2' |
1559 |
|
|
! ENDDO |
1560 |
|
|
! ENDIF |
1561 |
|
|
|
1562 |
|
|
|
1563 |
|
|
!<LTP |
1564 |
|
|
|
1565 |
|
|
!Limitation de la fraction surfacique couverte par les précipitations lorsque l'intensité locale du flux de précipitation descend en |
1566 |
|
|
!dessous de rain_int_min |
1567 |
|
|
IF (iflag_evap_prec==4) THEN |
1568 |
|
|
DO i=1, klon |
1569 |
|
|
IF (zrflclr(i) + ziflclr(i) .GT. 0 ) THEN |
1570 |
|
|
znebprecipclr(i) = min(znebprecipclr(i), max(zrflclr(i) & |
1571 |
|
|
/(znebprecipclr(i)*rain_int_min), & |
1572 |
|
|
ziflclr(i)/(znebprecipclr(i)*rain_int_min))) |
1573 |
|
|
ELSE |
1574 |
|
|
znebprecipclr(i)=0. |
1575 |
|
|
ENDIF |
1576 |
|
|
|
1577 |
|
|
IF (zrflcld(i) + ziflcld(i) .GT. 0 ) THEN |
1578 |
|
|
znebprecipcld(i) = min(znebprecipcld(i), & |
1579 |
|
|
max(zrflcld(i)/(znebprecipcld(i)*rain_int_min), & |
1580 |
|
|
ziflcld(i)/(znebprecipcld(i)*rain_int_min))) |
1581 |
|
|
ELSE |
1582 |
|
|
znebprecipcld(i)=0. |
1583 |
|
|
ENDIF |
1584 |
|
|
ENDDO |
1585 |
|
|
ENDIf |
1586 |
|
|
|
1587 |
|
|
!>LTP |
1588 |
|
|
|
1589 |
|
|
|
1590 |
|
|
|
1591 |
|
|
|
1592 |
|
|
|
1593 |
|
|
IF (.NOT. ice_thermo) THEN |
1594 |
|
|
DO i = 1, klon |
1595 |
|
|
IF (zt(i).LT.RTT) THEN |
1596 |
|
|
psfl(i,k)=zrfl(i) |
1597 |
|
|
ELSE |
1598 |
|
|
prfl(i,k)=zrfl(i) |
1599 |
|
|
ENDIF |
1600 |
|
|
ENDDO |
1601 |
|
|
ELSE |
1602 |
|
|
! JAM************************************************* |
1603 |
|
|
! Revoir partie ci-dessous: a quoi servent psfl et prfl? |
1604 |
|
|
! ***************************************************** |
1605 |
|
|
|
1606 |
|
|
DO i = 1, klon |
1607 |
|
|
! IF (zt(i).LT.RTT) THEN |
1608 |
|
|
psfl(i,k)=zifl(i) |
1609 |
|
|
! ELSE |
1610 |
|
|
prfl(i,k)=zrfl(i) |
1611 |
|
|
! ENDIF |
1612 |
|
|
!>AJ |
1613 |
|
|
ENDDO |
1614 |
|
|
ENDIF |
1615 |
|
|
! ---------------------------------------------------------------- |
1616 |
|
|
! Fin de formation des precipitations |
1617 |
|
|
! ---------------------------------------------------------------- |
1618 |
|
|
! |
1619 |
|
|
! Calculer les tendances de q et de t: |
1620 |
|
|
! |
1621 |
|
|
DO i = 1, klon |
1622 |
|
|
d_q(i,k) = zq(i) - q(i,k) |
1623 |
|
|
d_t(i,k) = zt(i) - t(i,k) |
1624 |
|
|
ENDDO |
1625 |
|
|
! |
1626 |
|
|
!AA--------------- Calcul du lessivage stratiforme ------------- |
1627 |
|
|
|
1628 |
|
|
DO i = 1,klon |
1629 |
|
|
! |
1630 |
|
|
if(zcond(i).gt.zoliq(i)+1.e-10) then |
1631 |
|
|
beta(i,k) = (zcond(i)-zoliq(i))/zcond(i)/dtime |
1632 |
|
|
else |
1633 |
|
|
beta(i,k) = 0. |
1634 |
|
|
endif |
1635 |
|
|
zprec_cond(i) = MAX(zcond(i)-zoliq(i),0.0) & |
1636 |
|
|
* (paprs(i,k)-paprs(i,k+1))/RG |
1637 |
|
|
IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN |
1638 |
|
|
!AA lessivage nucleation LMD5 dans la couche elle-meme |
1639 |
|
|
IF (iflag_t_glace.EQ.0) THEN |
1640 |
|
|
if (t(i,k) .GE. t_glace_min_old) THEN |
1641 |
|
|
zalpha_tr = a_tr_sca(3) |
1642 |
|
|
else |
1643 |
|
|
zalpha_tr = a_tr_sca(4) |
1644 |
|
|
endif |
1645 |
|
|
ELSE ! of IF (iflag_t_glace.EQ.0) |
1646 |
|
|
if (t(i,k) .GE. t_glace_min) THEN |
1647 |
|
|
zalpha_tr = a_tr_sca(3) |
1648 |
|
|
else |
1649 |
|
|
zalpha_tr = a_tr_sca(4) |
1650 |
|
|
endif |
1651 |
|
|
ENDIF |
1652 |
|
|
zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i)) |
1653 |
|
|
pfrac_nucl(i,k)=pfrac_nucl(i,k)*(1.-zneb(i)*zfrac_lessi) |
1654 |
|
|
frac_nucl(i,k)= 1.-zneb(i)*zfrac_lessi |
1655 |
|
|
! |
1656 |
|
|
! nucleation avec un facteur -1 au lieu de -0.5 |
1657 |
|
|
zfrac_lessi = 1. - EXP(-zprec_cond(i)/zneb(i)) |
1658 |
|
|
pfrac_1nucl(i,k)=pfrac_1nucl(i,k)*(1.-zneb(i)*zfrac_lessi) |
1659 |
|
|
ENDIF |
1660 |
|
|
! |
1661 |
|
|
ENDDO ! boucle sur i |
1662 |
|
|
! |
1663 |
|
|
!AA Lessivage par impaction dans les couches en-dessous |
1664 |
|
|
DO kk = k-1, 1, -1 |
1665 |
|
|
DO i = 1, klon |
1666 |
|
|
IF (rneb(i,k).GT.0.0.and.zprec_cond(i).gt.0.) THEN |
1667 |
|
|
IF (iflag_t_glace.EQ.0) THEN |
1668 |
|
|
if (t(i,kk) .GE. t_glace_min_old) THEN |
1669 |
|
|
zalpha_tr = a_tr_sca(1) |
1670 |
|
|
else |
1671 |
|
|
zalpha_tr = a_tr_sca(2) |
1672 |
|
|
endif |
1673 |
|
|
ELSE ! of IF (iflag_t_glace.EQ.0) |
1674 |
|
|
if (t(i,kk) .GE. t_glace_min) THEN |
1675 |
|
|
zalpha_tr = a_tr_sca(1) |
1676 |
|
|
else |
1677 |
|
|
zalpha_tr = a_tr_sca(2) |
1678 |
|
|
endif |
1679 |
|
|
ENDIF |
1680 |
|
|
zfrac_lessi = 1. - EXP(zalpha_tr*zprec_cond(i)/zneb(i)) |
1681 |
|
|
pfrac_impa(i,kk)=pfrac_impa(i,kk)*(1.-zneb(i)*zfrac_lessi) |
1682 |
|
|
frac_impa(i,kk)= 1.-zneb(i)*zfrac_lessi |
1683 |
|
|
ENDIF |
1684 |
|
|
ENDDO |
1685 |
|
|
ENDDO |
1686 |
|
|
! |
1687 |
|
|
!AA=============================================================== |
1688 |
|
|
! FIN DE LA BOUCLE VERTICALE |
1689 |
|
|
end DO |
1690 |
|
|
! |
1691 |
|
|
!AA================================================================== |
1692 |
|
|
! |
1693 |
|
|
! Pluie ou neige au sol selon la temperature de la 1ere couche |
1694 |
|
|
! |
1695 |
|
|
!CR: si la thermo de la glace est active, on calcule zifl directement |
1696 |
|
|
IF (.NOT.ice_thermo) THEN |
1697 |
|
|
DO i = 1, klon |
1698 |
|
|
IF ((t(i,1)+d_t(i,1)) .LT. RTT) THEN |
1699 |
|
|
!AJ< |
1700 |
|
|
! snow(i) = zrfl(i) |
1701 |
|
|
snow(i) = zrfl(i)+zifl(i) |
1702 |
|
|
!>AJ |
1703 |
|
|
zlh_solid(i) = RLSTT-RLVTT |
1704 |
|
|
ELSE |
1705 |
|
|
rain(i) = zrfl(i) |
1706 |
|
|
zlh_solid(i) = 0. |
1707 |
|
|
ENDIF |
1708 |
|
|
ENDDO |
1709 |
|
|
|
1710 |
|
|
ELSE |
1711 |
|
|
DO i = 1, klon |
1712 |
|
|
snow(i) = zifl(i) |
1713 |
|
|
rain(i) = zrfl(i) |
1714 |
|
|
ENDDO |
1715 |
|
|
|
1716 |
|
|
ENDIF |
1717 |
|
|
! |
1718 |
|
|
! For energy conservation : when snow is present, the solification |
1719 |
|
|
! latent heat is considered. |
1720 |
|
|
!CR: si thermo de la glace, neige deja prise en compte |
1721 |
|
|
IF (.not.ice_thermo) THEN |
1722 |
|
|
DO k = 1, klev |
1723 |
|
|
DO i = 1, klon |
1724 |
|
|
zcpair=RCPD*(1.0+RVTMP2*(q(i,k)+d_q(i,k))) |
1725 |
|
|
zmair(i)=(paprs(i,k)-paprs(i,k+1))/RG |
1726 |
|
|
zm_solid = (prfl(i,k)-prfl(i,k+1)+psfl(i,k)-psfl(i,k+1))*dtime |
1727 |
|
|
d_t(i,k) = d_t(i,k) + zlh_solid(i) *zm_solid / (zcpair*zmair(i)) |
1728 |
|
|
END DO |
1729 |
|
|
END DO |
1730 |
|
|
ENDIF |
1731 |
|
|
! |
1732 |
|
|
|
1733 |
|
|
if (ncoreczq>0) then |
1734 |
|
|
WRITE(lunout,*)'WARNING : ZQ dans fisrtilp ',ncoreczq,' val < 1.e-15.' |
1735 |
|
|
endif |
1736 |
|
|
|
1737 |
|
|
END SUBROUTINE fisrtilp |