1 |
|
|
! |
2 |
|
|
! $Id $ |
3 |
|
|
! |
4 |
|
|
SUBROUTINE cvltr_spl(pdtime, da, phi,phi2,d1a,dam, mpIN,epIN, & |
5 |
|
|
sigd,sij,wght_cvfd,clw,elij,epmlmMm,eplaMm, & |
6 |
|
|
pmflxrIN,pmflxsIN,ev,te,wdtrainA,wdtrainM, & |
7 |
|
|
paprs,it,tr,upd,dnd,inb,icb, & |
8 |
|
|
kk,henry,zrho, ccntrAA_spla,ccntrENV_spla,coefcoli_spla, & |
9 |
|
|
id_prec,id_fine,id_coss, id_codu, id_scdu, & |
10 |
|
|
dtrcv,trsptd,dtrSscav,dtrsat,dtrUscav,qDi,qPr, & |
11 |
|
|
qPa,qMel,qTrdi,dtrcvMA,Mint, & |
12 |
|
|
zmfd1a,zmfphi2,zmfdam) |
13 |
|
|
USE IOIPSL |
14 |
|
|
USE dimphy |
15 |
|
|
USE infotrac_phy, ONLY : nbtr |
16 |
|
|
IMPLICIT NONE |
17 |
|
|
!===================================================================== |
18 |
|
|
! Objet : convection des traceurs / KE |
19 |
|
|
! Auteurs: M-A Filiberti and J-Y Grandpeix |
20 |
|
|
! modifiee par R Pilon : lessivage des traceurs / KE |
21 |
|
|
!===================================================================== |
22 |
|
|
|
23 |
|
|
include "YOMCST.h" |
24 |
|
|
include "YOECUMF.h" |
25 |
|
|
include "conema3.h" |
26 |
|
|
include "chem.h" |
27 |
|
|
|
28 |
|
|
! Entree |
29 |
|
|
REAL,INTENT(IN) :: pdtime |
30 |
|
|
REAL,DIMENSION(klon,klev),INTENT(IN) :: da |
31 |
|
|
REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi |
32 |
|
|
! RomP |
33 |
|
|
REAL,DIMENSION(klon,klev),INTENT(IN) :: d1a,dam ! matrices pour simplifier |
34 |
|
|
REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi2 ! l'ecriture des tendances |
35 |
|
|
! |
36 |
|
|
REAL,DIMENSION(klon,klev),INTENT(IN) :: mpIN |
37 |
|
|
REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression aux 1/2 couches (bas en haut) |
38 |
|
|
! REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression aux 1/2 couches (bas en haut) |
39 |
|
|
REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: tr ! q de traceur (bas en haut) |
40 |
|
|
INTEGER,INTENT(IN) :: it |
41 |
|
|
REAL,DIMENSION(klon,klev),INTENT(IN) :: upd ! saturated updraft mass flux |
42 |
|
|
REAL,DIMENSION(klon,klev),INTENT(IN) :: dnd ! saturated downdraft mass flux |
43 |
|
|
! |
44 |
|
|
REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtrainA ! masses precipitantes de l'asc adiab |
45 |
|
|
REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtrainM ! masses precipitantes des melanges |
46 |
|
|
!JE REAL,DIMENSION(klon,klev),INTENT(IN) :: pmflxrIN ! vprecip: eau |
47 |
|
|
REAL,DIMENSION(klon,klev+1),INTENT(IN) :: pmflxrIN ! vprecip: eau |
48 |
|
|
!JE REAL,DIMENSION(klon,klev),INTENT(IN) :: pmflxsIN ! vprecip: neige |
49 |
|
|
REAL,DIMENSION(klon,klev+1),INTENT(IN) :: pmflxsIN ! vprecip: neige |
50 |
|
|
REAL,DIMENSION(klon,klev),INTENT(IN) :: ev ! evaporation cv30_routine |
51 |
|
|
REAL,DIMENSION(klon,klev),INTENT(IN) :: epIN |
52 |
|
|
REAL,DIMENSION(klon,klev),INTENT(IN) :: te |
53 |
|
|
REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: sij ! fraction dair de lenv |
54 |
|
|
REAL,DIMENSION(klon,klev),INTENT(IN) :: wght_cvfd ! weights of the layers feeding convection |
55 |
|
|
REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: elij ! contenu en eau condensée spécifique/conc deau condensée massique |
56 |
|
|
REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: epmlmMm ! eau condensee precipitee dans mel masse dair sat |
57 |
|
|
REAL,DIMENSION(klon,klev),INTENT(IN) :: eplaMm ! eau condensee precipitee dans aa masse dair sat |
58 |
|
|
|
59 |
|
|
REAL,DIMENSION(klon,klev),INTENT(IN) :: clw ! contenu en eau condensée dans lasc adiab |
60 |
|
|
REAL,DIMENSION(klon),INTENT(IN) :: sigd |
61 |
|
|
INTEGER,DIMENSION(klon),INTENT(IN) :: icb,inb |
62 |
|
|
! Sortie |
63 |
|
|
REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrcv ! tendance totale (bas en haut) |
64 |
|
|
REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrcvMA ! M-A Filiberti |
65 |
|
|
REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: trsptd ! tendance du transport |
66 |
|
|
REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrSscav ! tendance du lessivage courant sat |
67 |
|
|
REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrsat ! tendance trsp+sat scav |
68 |
|
|
REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: dtrUscav ! tendance du lessivage courant unsat |
69 |
|
|
! |
70 |
|
|
! Variables locales |
71 |
|
|
INTEGER :: i,j,k |
72 |
|
|
REAL,DIMENSION(klon,klev) :: dxpres ! difference de pression entre niveau (j+1) et (j) |
73 |
|
|
REAL :: pdtimeRG ! pas de temps * gravite |
74 |
|
|
REAL,DIMENSION(klon,nbtr) :: qfeed ! tracer concentration feeding convection |
75 |
|
|
! variables pour les courants satures |
76 |
|
|
REAL,DIMENSION(klon,klev,klev) :: zmd |
77 |
|
|
REAL,DIMENSION(klon,klev,klev) :: za |
78 |
|
|
REAL,DIMENSION(klon,klev,nbtr) :: zmfd,zmfa |
79 |
|
|
REAL,DIMENSION(klon,klev,nbtr) :: zmfp,zmfu |
80 |
|
|
|
81 |
|
|
REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfd1a |
82 |
|
|
REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfdam |
83 |
|
|
REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: zmfphi2 |
84 |
|
|
|
85 |
|
|
! RomP ! les variables sont nettoyees des valeurs aberrantes |
86 |
|
|
REAL,DIMENSION(klon,klev) :: Pa, Pm ! pluie AA et mélanges, var temporaire |
87 |
|
|
REAL,DIMENSION(klon,klev) :: pmflxs,pmflxr ! pmflxrIN,pmflxsIN sans valeur aberante |
88 |
|
|
REAL,DIMENSION(klon,klev) :: mp ! flux de masse |
89 |
|
|
REAL,DIMENSION(klon,klev) :: ep ! fraction d'eau convertie en precipitation |
90 |
|
|
REAL,DIMENSION(klon,klev) :: evap ! evaporation : variable temporaire |
91 |
|
|
REAL,DIMENSION(klon,klev) :: rho !environmental density |
92 |
|
|
|
93 |
|
|
REAL,DIMENSION(klon,klev) :: kappa ! denominateur du au calcul de la matrice |
94 |
|
|
! pour obtenir qd et qp |
95 |
|
|
REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qTrdi ! traceurs descente air insature transport MA |
96 |
|
|
REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qDi ! traceurs descente insaturees |
97 |
|
|
REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qPr ! traceurs colonne precipitante |
98 |
|
|
REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qPa ! traceurs dans les precip issues lasc. adiab. |
99 |
|
|
REAL,DIMENSION(klon,klev,nbtr),INTENT(OUT) :: qMel ! traceurs dans les precip issues des melanges |
100 |
|
|
REAL,DIMENSION(klon,klev,nbtr) :: qMeltmp ! variable temporaire |
101 |
|
|
REAL,DIMENSION(klon,klev,nbtr) :: qpmMint |
102 |
|
|
REAL,DIMENSION(klon,klev),INTENT(OUT) :: Mint |
103 |
|
|
! tendances |
104 |
|
|
REAL :: tdcvMA ! terme de transport de traceur (schema Marie Angele) |
105 |
|
|
REAL :: trsptrac ! terme de transport de traceur par l'air |
106 |
|
|
REAL :: scavtrac ! terme de lessivage courant sature |
107 |
|
|
REAL :: uscavtrac ! terme de lessivage courant insature |
108 |
|
|
! impaction |
109 |
|
|
!!! Correction apres discussion Romain P. / Olivier B. |
110 |
|
|
!!! REAL,PARAMETER :: rdrop=2.5e-3 ! rayon des gouttes d'eau |
111 |
|
|
REAL,PARAMETER :: rdrop=1.e-3 ! rayon des gouttes d'eau |
112 |
|
|
!!! |
113 |
|
|
REAL,DIMENSION(klon,klev) :: imp ! coefficient d'impaction |
114 |
|
|
! parametres lessivage |
115 |
|
|
REAL :: ccntrAA_coef ! \alpha_a : fract aerosols de l'AA convertis en CCN |
116 |
|
|
REAL :: ccntrENV_coef ! \beta_m : fract aerosols de l'env convertis en CCN |
117 |
|
|
REAL :: coefcoli ! coefficient de collision des gouttes sur les aerosols |
118 |
|
|
! |
119 |
|
|
LOGICAL,DIMENSION(klon,klev) :: NO_precip |
120 |
|
|
! LOGICAL :: scavON |
121 |
|
|
! var tmp tests |
122 |
|
|
REAL :: conserv |
123 |
|
|
real :: conservMA |
124 |
|
|
! JE SPLA adaptation |
125 |
|
|
|
126 |
|
|
INTEGER :: id_prec,id_fine,id_coss, id_codu, id_scdu |
127 |
|
|
REAL,DIMENSION(nbtr) :: ccntrAA_spla,ccntrENV_spla,coefcoli_spla |
128 |
|
|
REAL,DIMENSION(klon,klev) :: ccntrAA_coef3d |
129 |
|
|
REAL,DIMENSION(klon,klev) :: ccntrENV_coef3d |
130 |
|
|
REAL,DIMENSION(klon,klev) :: coefcoli3d |
131 |
|
|
|
132 |
|
|
REAL,DIMENSION(nbtr) :: henry !--cste de Henry mol/l/atm |
133 |
|
|
REAL,DIMENSION(nbtr) :: kk !--coefficient de var avec T (K) |
134 |
|
|
REAL :: henry_t !--constante de Henry a T t (mol/l/atm) |
135 |
|
|
REAL :: f_a !--rapport de la phase aqueuse a la phase gazeuse |
136 |
|
|
|
137 |
|
|
REAL, PARAMETER :: ph=5. |
138 |
|
|
REAL :: K1, K2 |
139 |
|
|
REAL,DIMENSION(klon,klev) :: zrho |
140 |
|
|
REAL, PARAMETER :: qliq=1.e-3 ! CONVECTIVE ONLY |
141 |
|
|
|
142 |
|
|
! Je end 20140105 |
143 |
|
|
|
144 |
|
|
!JE20140724<< |
145 |
|
|
!JE SPLA adapt: |
146 |
|
|
! |
147 |
|
|
! |
148 |
|
|
!! coefficient lessivage |
149 |
|
|
! ccntrAA_coef = 0. |
150 |
|
|
! ccntrENV_coef = 0. |
151 |
|
|
! coefcoli = 0. |
152 |
|
|
! |
153 |
|
|
!!$OMP MASTER |
154 |
|
|
! call getin('ccntrAA_coef',ccntrAA_coef) |
155 |
|
|
! call getin('ccntrENV_coef',ccntrENV_coef) |
156 |
|
|
! call getin('coefcoli',coefcoli) |
157 |
|
|
!!$OMP END MASTER |
158 |
|
|
!!$OMP BARRIER |
159 |
|
|
!!! JE |
160 |
|
|
! do j=1,klev |
161 |
|
|
! do i=1,klon |
162 |
|
|
! imp(i,j)=0. |
163 |
|
|
! ccntrAA_coef3d(i,j)=ccntrAA_coef |
164 |
|
|
! ccntrENV_coef3d(i,j)=ccntrENV_coef |
165 |
|
|
! coefcoli3d(i,j)=coefcoli |
166 |
|
|
! enddo |
167 |
|
|
! enddo |
168 |
|
|
! |
169 |
|
|
! |
170 |
|
|
!!! for SPLA |
171 |
|
|
!!!JEend |
172 |
|
|
! print*,'cvltr coef lessivage convectif', ccntrAA_coef,ccntrENV_coef,coefcoli |
173 |
|
|
! |
174 |
|
|
!JE20140724>> |
175 |
|
|
|
176 |
|
|
! scavON=.TRUE. |
177 |
|
|
! if(scavON) then |
178 |
|
|
! ccntrAA_coef = 1. |
179 |
|
|
! ccntrENV_coef = 1. |
180 |
|
|
! coefcoli = 1. |
181 |
|
|
! else |
182 |
|
|
! ccntrAA_coef = 0. |
183 |
|
|
! ccntrENV_coef = 0. |
184 |
|
|
! coefcoli = 0. |
185 |
|
|
! endif |
186 |
|
|
|
187 |
|
|
! ====================================================== |
188 |
|
|
! calcul de l'impaction |
189 |
|
|
! ====================================================== |
190 |
|
|
!initialisation |
191 |
|
|
do j=1,klev |
192 |
|
|
do i=1,klon |
193 |
|
|
imp(i,j)=0. |
194 |
|
|
enddo |
195 |
|
|
enddo |
196 |
|
|
!JE init 20140103 |
197 |
|
|
!! impaction sur la surface de la colonne de la descente insaturee |
198 |
|
|
!! On prend la moyenne des precip entre le niveau i+1 et i |
199 |
|
|
!! I=3/4* (P(1+1)+P(i))/2 / (sigd*r*rho_l) |
200 |
|
|
!! 1000kg/m3= densité de l'eau |
201 |
|
|
!! 0.75e-3 = 3/4 /1000 |
202 |
|
|
!! Par la suite, I est tout le temps multiplié par sig_d pour avoir l'impaction sur la surface de la maille |
203 |
|
|
!! on le néglige ici pour simplifier le code |
204 |
|
|
! do j=1,klev-1 |
205 |
|
|
! do i=1,klon |
206 |
|
|
! imp(i,j) = coefcoli*0.75e-3/rdrop *& |
207 |
|
|
! 0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j)) |
208 |
|
|
!! rho(i,j)=pplay(i,j)/(rd*te(i,j)) |
209 |
|
|
! enddo |
210 |
|
|
! enddo |
211 |
|
|
!JE end 20140103 |
212 |
|
|
|
213 |
|
|
! |
214 |
|
|
! initialisation pour flux de traceurs, td et autre |
215 |
|
|
trsptrac = 0. |
216 |
|
|
scavtrac = 0. |
217 |
|
|
uscavtrac = 0. |
218 |
|
|
qfeed(:,it) = 0. !RL |
219 |
|
|
DO j=1,klev |
220 |
|
|
DO i=1,klon |
221 |
|
|
zmfd(i,j,it)=0. |
222 |
|
|
zmfa(i,j,it)=0. |
223 |
|
|
zmfu(i,j,it)=0. |
224 |
|
|
zmfp(i,j,it)=0. |
225 |
|
|
zmfphi2(i,j,it)=0. |
226 |
|
|
zmfd1a(i,j,it)=0. |
227 |
|
|
zmfdam(i,j,it)=0. |
228 |
|
|
qDi(i,j,it)=0. |
229 |
|
|
qPr(i,j,it)=0. |
230 |
|
|
qPa(i,j,it)=0. |
231 |
|
|
qMel(i,j,it)=0. |
232 |
|
|
qMeltmp(i,j,it)=0. |
233 |
|
|
qTrdi(i,j,it)=0. |
234 |
|
|
kappa(i,j)=0. |
235 |
|
|
trsptd(i,j,it)=0. |
236 |
|
|
dtrsat(i,j,it)=0. |
237 |
|
|
dtrSscav(i,j,it)=0. |
238 |
|
|
dtrUscav(i,j,it)=0. |
239 |
|
|
dtrcv(i,j,it)=0. |
240 |
|
|
dtrcvMA(i,j,it)=0. |
241 |
|
|
evap(i,j)=0. |
242 |
|
|
dxpres(i,j)=0. |
243 |
|
|
qpmMint(i,j,it)=0. |
244 |
|
|
Mint(i,j)=0. |
245 |
|
|
END DO |
246 |
|
|
END DO |
247 |
|
|
|
248 |
|
|
! suppression des valeurs très faibles (~1e-320) |
249 |
|
|
! multiplication de levaporation pour lavoir par unite de temps |
250 |
|
|
! et par unite de surface de la maille |
251 |
|
|
! -> cv30_unsat : evap : masse evaporee/s/(m2 de la descente) |
252 |
|
|
DO j=1,klev |
253 |
|
|
DO i=1,klon |
254 |
|
|
if(ev(i,j).lt.1.e-16) then |
255 |
|
|
evap(i,j)=0. |
256 |
|
|
else |
257 |
|
|
evap(i,j)=ev(i,j)*sigd(i) |
258 |
|
|
endif |
259 |
|
|
END DO |
260 |
|
|
END DO |
261 |
|
|
|
262 |
|
|
DO j=1,klev |
263 |
|
|
DO i=1,klon |
264 |
|
|
if(j.lt.klev) then |
265 |
|
|
if(epIN(i,j).lt.1.e-32) then |
266 |
|
|
ep(i,j)=0. |
267 |
|
|
else |
268 |
|
|
ep(i,j)=epIN(i,j) |
269 |
|
|
endif |
270 |
|
|
else |
271 |
|
|
ep(i,j)=epmax |
272 |
|
|
endif |
273 |
|
|
if(mpIN(i,j).lt.1.e-32) then |
274 |
|
|
mp(i,j)=0. |
275 |
|
|
else |
276 |
|
|
mp(i,j)=mpIN(i,j) |
277 |
|
|
endif |
278 |
|
|
if(pmflxsIN(i,j).lt.1.e-32) then |
279 |
|
|
pmflxs(i,j)=0. |
280 |
|
|
else |
281 |
|
|
pmflxs(i,j)=pmflxsIN(i,j) |
282 |
|
|
endif |
283 |
|
|
if(pmflxrIN(i,j).lt.1.e-32) then |
284 |
|
|
pmflxr(i,j)=0. |
285 |
|
|
else |
286 |
|
|
pmflxr(i,j)=pmflxrIN(i,j) |
287 |
|
|
endif |
288 |
|
|
if(wdtrainA(i,j).lt.1.e-32) then |
289 |
|
|
Pa(i,j)=0. |
290 |
|
|
else |
291 |
|
|
Pa(i,j)=wdtrainA(i,j) |
292 |
|
|
endif |
293 |
|
|
if(wdtrainM(i,j).lt.1.e-32) then |
294 |
|
|
Pm(i,j)=0. |
295 |
|
|
else |
296 |
|
|
Pm(i,j)=wdtrainM(i,j) |
297 |
|
|
endif |
298 |
|
|
END DO |
299 |
|
|
END DO |
300 |
|
|
|
301 |
|
|
!========================================== |
302 |
|
|
DO j = klev-1,1,-1 |
303 |
|
|
DO i = 1,klon |
304 |
|
|
NO_precip(i,j) = (pmflxr(i,j+1)+pmflxs(i,j+1)).lt.1.e-10& |
305 |
|
|
.and.Pa(i,j).lt.1.e-10.and.Pm(i,j).lt.1.e-10 |
306 |
|
|
END DO |
307 |
|
|
END DO |
308 |
|
|
!============================================================================== |
309 |
|
|
! JE SPLA: Calc of ccntrAA_coef,ccntrENV_coef, coefcoli for SPLA aerosols and |
310 |
|
|
! precursors. From SPLA code inscav_spl.F |
311 |
|
|
!print *,'Using SPLA values for cvltr_spl ccntr and coefcoli params' |
312 |
|
|
! |
313 |
|
|
! |
314 |
|
|
!IF (it.EQ.2) THEN !--fine aerosol |
315 |
|
|
! DO j=1,klev |
316 |
|
|
! DO i=1,klon |
317 |
|
|
! ccntrAA_coef3d(i,j)=0.7 |
318 |
|
|
! ccntrENV_coef3d(i,j)=0.7 |
319 |
|
|
! coefcoli3d(i,j)=0.001 |
320 |
|
|
! ENDDO |
321 |
|
|
! ENDDO |
322 |
|
|
!ELSEIF (it.EQ.3) THEN !-- Coarse Sea salt aerosol |
323 |
|
|
! DO j=1,klev |
324 |
|
|
! DO i=1,klon |
325 |
|
|
! ccntrAA_coef3d(i,j)=1.0 |
326 |
|
|
! ccntrENV_coef3d(i,j)=1.0 |
327 |
|
|
! coefcoli3d(i,j)=0.001 |
328 |
|
|
! ENDDO |
329 |
|
|
! ENDDO |
330 |
|
|
! |
331 |
|
|
!ELSEIF (it.EQ.4) THEN !--Coarse Dust aerosol |
332 |
|
|
! DO j=1,klev |
333 |
|
|
! DO i=1,klon |
334 |
|
|
! ccntrAA_coef3d(i,j)=0.7 |
335 |
|
|
! ccntrENV_coef3d(i,j)=0.7 |
336 |
|
|
! coefcoli3d(i,j)=0.001 |
337 |
|
|
! |
338 |
|
|
! ENDDO |
339 |
|
|
! ENDDO |
340 |
|
|
! Gas precursor: Henry's law |
341 |
|
|
|
342 |
|
|
IF (it .EQ. id_prec) THEN |
343 |
|
|
DO k=1, klev |
344 |
|
|
DO i=1, klon |
345 |
|
|
henry_t=henry(it)*exp(-kk(it)*(1./298.-1./te(i,k))) !--mol/l/atm |
346 |
|
|
K1=1.2e-2*exp(-2010*(1/298.-1/te(i,k))) |
347 |
|
|
K2=6.6e-8*exp(-1510*(1/298.-1/te(i,k))) |
348 |
|
|
henry_t=henry_t*(1. + K1/10.**(-ph) + K1*K2/(10.**(-ph))**2) |
349 |
|
|
f_a=henry_t/101.325*R*te(i,k)*qliq*zrho(i,k)/rho_water |
350 |
|
|
! scav(i,k)=f_a/(1.+f_a) |
351 |
|
|
ccntrAA_coef3d(i,k)= f_a/(1.+f_a) |
352 |
|
|
ccntrENV_coef3d(i,k)= f_a/(1.+f_a) |
353 |
|
|
coefcoli3d(i,k)=0.0 |
354 |
|
|
ENDDO |
355 |
|
|
ENDDO |
356 |
|
|
! CALL minmaxqfi2(clw,1.e33,-1.e33,'minmax clw') |
357 |
|
|
ELSE |
358 |
|
|
DO j=1,klev |
359 |
|
|
DO i=1,klon |
360 |
|
|
ccntrAA_coef3d(i,j)=ccntrAA_spla(it) |
361 |
|
|
ccntrENV_coef3d(i,j)=ccntrENV_spla(it) |
362 |
|
|
coefcoli3d(i,j)=coefcoli_spla(it) |
363 |
|
|
ENDDO |
364 |
|
|
ENDDO |
365 |
|
|
|
366 |
|
|
|
367 |
|
|
ENDIF |
368 |
|
|
|
369 |
|
|
! JE end SPLA modifs in ccn parameters |
370 |
|
|
!============================================================================== |
371 |
|
|
|
372 |
|
|
|
373 |
|
|
|
374 |
|
|
|
375 |
|
|
!JE init 20140103 |
376 |
|
|
! impaction sur la surface de la colonne de la descente insaturee |
377 |
|
|
! On prend la moyenne des precip entre le niveau i+1 et i |
378 |
|
|
! I=3/4* (P(1+1)+P(i))/2 / (sigd*r*rho_l) |
379 |
|
|
! 1000kg/m3= densité de l'eau |
380 |
|
|
! 0.75e-3 = 3/4 /1000 |
381 |
|
|
! Par la suite, I est tout le temps multiplié par sig_d pour avoir l'impaction sur la surface de la maille |
382 |
|
|
! on le néglige ici pour simplifier le code |
383 |
|
|
do j=1,klev-1 |
384 |
|
|
do i=1,klon |
385 |
|
|
!JE imp(i,j) = coefcoli*0.75e-3/rdrop *& |
386 |
|
|
!JE 0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j)) |
387 |
|
|
imp(i,j) = coefcoli3d(i,j)*0.75e-3/rdrop *& |
388 |
|
|
0.5*(pmflxr(i,j+1)+pmflxs(i,j+1)+pmflxr(i,j)+pmflxs(i,j)) |
389 |
|
|
! rho(i,j)=pplay(i,j)/(rd*te(i,j)) |
390 |
|
|
enddo |
391 |
|
|
enddo |
392 |
|
|
!JE end 20140103 |
393 |
|
|
|
394 |
|
|
|
395 |
|
|
! ========================================= |
396 |
|
|
! calcul des tendances liees au downdraft |
397 |
|
|
! ========================================= |
398 |
|
|
!cdir collapse |
399 |
|
|
DO k=1,klev |
400 |
|
|
DO j=1,klev |
401 |
|
|
DO i=1,klon |
402 |
|
|
zmd(i,j,k)=0. |
403 |
|
|
za (i,j,k)=0. |
404 |
|
|
END DO |
405 |
|
|
END DO |
406 |
|
|
END DO |
407 |
|
|
! calcul de la matrice d echange |
408 |
|
|
! matrice de distribution de la masse entrainee en k |
409 |
|
|
! commmentaire RomP : mp > 0 |
410 |
|
|
DO k=1,klev-1 |
411 |
|
|
DO i=1,klon |
412 |
|
|
zmd(i,k,k)=max(0.,mp(i,k)-mp(i,k+1)) ! ~ mk(k) |
413 |
|
|
END DO |
414 |
|
|
END DO |
415 |
|
|
DO k=2,klev |
416 |
|
|
DO j=k-1,1,-1 |
417 |
|
|
DO i=1,klon |
418 |
|
|
if(mp(i,j+1).gt.1.e-10) then |
419 |
|
|
zmd(i,j,k)=zmd(i,j+1,k)*min(1.,mp(i,j)/mp(i,j+1)) !det ~ mk(j)=mk(j+1)*mp(i,j)/mp(i,j+1) |
420 |
|
|
ENDif |
421 |
|
|
END DO |
422 |
|
|
END DO |
423 |
|
|
END DO |
424 |
|
|
DO k=1,klev |
425 |
|
|
DO j=1,klev-1 |
426 |
|
|
DO i=1,klon |
427 |
|
|
za(i,j,k)=max(0.,zmd(i,j+1,k)-zmd(i,j,k)) |
428 |
|
|
END DO |
429 |
|
|
END DO |
430 |
|
|
END DO |
431 |
|
|
!!!!! quantite de traceur dans la descente d'air insaturee : 4 juin 2012 |
432 |
|
|
DO k=1,klev |
433 |
|
|
DO j=1,klev-1 |
434 |
|
|
DO i=1,klon |
435 |
|
|
if(mp(i,j+1).gt.1.e-10) then |
436 |
|
|
qTrdi(i,j+1,it)=qTrdi(i,j+1,it)+(zmd(i,j+1,k)/mp(i,j+1))*tr(i,k,it) |
437 |
|
|
else |
438 |
|
|
qTrdi(i,j,it)=0.!tr(i,j,it) |
439 |
|
|
endif |
440 |
|
|
ENDDO |
441 |
|
|
ENDDO |
442 |
|
|
ENDDO |
443 |
|
|
!!!!! |
444 |
|
|
! |
445 |
|
|
! rajout du terme lie a l ascendance induite |
446 |
|
|
! |
447 |
|
|
DO j=2,klev |
448 |
|
|
DO i=1,klon |
449 |
|
|
za(i,j,j-1)=za(i,j,j-1)+mp(i,j) |
450 |
|
|
END DO |
451 |
|
|
END DO |
452 |
|
|
! |
453 |
|
|
! tendance courants insatures ! sans lessivage ancien schema |
454 |
|
|
! |
455 |
|
|
DO k=1,klev |
456 |
|
|
DO j=1,klev |
457 |
|
|
DO i=1,klon |
458 |
|
|
zmfd(i,j,it)=zmfd(i,j,it)+za(i,j,k)*(tr(i,k,it)-tr(i,j,it)) |
459 |
|
|
END DO |
460 |
|
|
END DO |
461 |
|
|
END DO |
462 |
|
|
! |
463 |
|
|
! ========================================= |
464 |
|
|
! calcul des tendances liees aux courants satures j <-> z ; k <-> z' |
465 |
|
|
! ========================================= |
466 |
|
|
!RL |
467 |
|
|
! Feeding concentrations |
468 |
|
|
DO j=1,klev |
469 |
|
|
DO i=1,klon |
470 |
|
|
qfeed(i,it)=qfeed(i,it)+wght_cvfd(i,j)*tr(i,j,it) |
471 |
|
|
END DO |
472 |
|
|
END DO |
473 |
|
|
!RL |
474 |
|
|
! |
475 |
|
|
DO j=1,klev |
476 |
|
|
DO i=1,klon |
477 |
|
|
!RL |
478 |
|
|
!! zmfa(i,j,it)=da(i,j)*(tr(i,1,it)-tr(i,j,it)) ! da |
479 |
|
|
zmfa(i,j,it)=da(i,j)*(qfeed(i,it)-tr(i,j,it)) ! da |
480 |
|
|
!RL |
481 |
|
|
END DO |
482 |
|
|
END DO |
483 |
|
|
! |
484 |
|
|
DO k=1,klev |
485 |
|
|
DO j=1,klev |
486 |
|
|
DO i=1,klon |
487 |
|
|
zmfp(i,j,it)=zmfp(i,j,it)+phi(i,j,k)*(tr(i,k,it)-tr(i,j,it)) ! phi |
488 |
|
|
END DO |
489 |
|
|
END DO |
490 |
|
|
END DO |
491 |
|
|
! RomP ajout des matrices liees au lessivage |
492 |
|
|
DO j=1,klev |
493 |
|
|
DO i=1,klon |
494 |
|
|
zmfd1a(i,j,it)=d1a(i,j)*tr(i,1,it) ! da1 |
495 |
|
|
zmfdam(i,j,it)=dam(i,j)*tr(i,1,it) ! dam |
496 |
|
|
END DO |
497 |
|
|
END DO |
498 |
|
|
DO k=1,klev |
499 |
|
|
DO j=1,klev |
500 |
|
|
DO i=1,klon |
501 |
|
|
zmfphi2(i,j,it)=zmfphi2(i,j,it)+phi2(i,j,k)*tr(i,k,it) ! psi |
502 |
|
|
END DO |
503 |
|
|
END DO |
504 |
|
|
END DO |
505 |
|
|
DO j=1,klev-1 |
506 |
|
|
DO i=1,klon |
507 |
|
|
zmfu(i,j,it)=max(0.,upd(i,j+1)+dnd(i,j+1))*(tr(i,j+1,it)-tr(i,j,it)) |
508 |
|
|
END DO |
509 |
|
|
END DO |
510 |
|
|
DO j=2,klev |
511 |
|
|
DO i=1,klon |
512 |
|
|
zmfu(i,j,it)=zmfu(i,j,it)+min(0.,upd(i,j)+dnd(i,j))*(tr(i,j,it)-tr(i,j-1,it)) |
513 |
|
|
END DO |
514 |
|
|
END DO |
515 |
|
|
! =================================================== |
516 |
|
|
! calcul des tendances liees aux courants insatures |
517 |
|
|
! =================================================== |
518 |
|
|
! pression |
519 |
|
|
DO k=1, klev |
520 |
|
|
DO i=1, klon |
521 |
|
|
dxpres(i,k)=paprs(i,k)-paprs(i,k+1) |
522 |
|
|
ENDDO |
523 |
|
|
ENDDO |
524 |
|
|
pdtimeRG=pdtime*RG |
525 |
|
|
|
526 |
|
|
! q_pa et q_pm traceurs issues des courants satures se retrouvant dans les precipitations |
527 |
|
|
DO j=1,klev |
528 |
|
|
DO i=1,klon |
529 |
|
|
if(j.ge.icb(i).and.j.le.inb(i)) then |
530 |
|
|
if(clw(i,j).gt.1.e-16) then |
531 |
|
|
! qPa(i,j,it)=ccntrAA_coef*tr(i,1,it)/clw(i,j) |
532 |
|
|
qPa(i,j,it)=ccntrAA_coef3d(i,j)*tr(i,1,it)/clw(i,j) |
533 |
|
|
else |
534 |
|
|
qPa(i,j,it)=0. |
535 |
|
|
endif |
536 |
|
|
endif |
537 |
|
|
END DO |
538 |
|
|
END DO |
539 |
|
|
|
540 |
|
|
! calcul de q_pm en 2 parties : |
541 |
|
|
! 1) calcul de sa valeur pour un niveau z' donne |
542 |
|
|
! 2) integration sur la verticale sur z' |
543 |
|
|
DO j=1,klev |
544 |
|
|
DO k=1,j-1 |
545 |
|
|
DO i=1,klon |
546 |
|
|
if(k.ge.icb(i).and.k.le.inb(i).and.& |
547 |
|
|
j.le.inb(i)) then |
548 |
|
|
if(elij(i,k,j).gt.1.e-16) then |
549 |
|
|
!JE qMeltmp(i,j,it)=((1-ep(i,k))*ccntrAA_coef*tr(i,1,it)& |
550 |
|
|
!JE *(1.-sij(i,k,j)) +ccntrENV_coef& |
551 |
|
|
!JE *tr(i,k,it)*sij(i,k,j)) / elij(i,k,j) |
552 |
|
|
qMeltmp(i,j,it)=((1-ep(i,k))*ccntrAA_coef3d(i,k)*tr(i,1,it)& |
553 |
|
|
*(1.-sij(i,k,j)) +ccntrENV_coef3d(i,k)& |
554 |
|
|
*tr(i,k,it)*sij(i,k,j)) / elij(i,k,j) |
555 |
|
|
else |
556 |
|
|
qMeltmp(i,j,it)=0. |
557 |
|
|
endif |
558 |
|
|
qpmMint(i,j,it)=qpmMint(i,j,it) + qMeltmp(i,j,it)*epmlmMm(i,j,k) |
559 |
|
|
Mint(i,j)=Mint(i,j) + epmlmMm(i,j,k) |
560 |
|
|
endif ! end if dans nuage |
561 |
|
|
END DO |
562 |
|
|
END DO |
563 |
|
|
END DO |
564 |
|
|
|
565 |
|
|
DO j=1,klev |
566 |
|
|
DO i=1,klon |
567 |
|
|
if(Mint(i,j).gt.1.e-16) then |
568 |
|
|
qMel(i,j,it)=qpmMint(i,j,it)/Mint(i,j) |
569 |
|
|
else |
570 |
|
|
qMel(i,j,it)=0. |
571 |
|
|
endif |
572 |
|
|
END DO |
573 |
|
|
END DO |
574 |
|
|
|
575 |
|
|
! calcul de q_d et q_p traceurs de la descente precipitante |
576 |
|
|
DO j=klev-1,1,-1 |
577 |
|
|
DO i=1,klon |
578 |
|
|
if(mp(i,j+1).gt.mp(i,j).and.mp(i,j+1).gt.1.e-10) then ! detrainement |
579 |
|
|
kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*& |
580 |
|
|
(-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))& |
581 |
|
|
+ (imp(i,j)/RG*dxpres(i,j))*(evap(i,j)/RG*dxpres(i,j))) |
582 |
|
|
|
583 |
|
|
elseif(mp(i,j).gt.mp(i,j+1).and.mp(i,j).gt.1.e-10) then! entrainement |
584 |
|
|
if(j.eq.1) then |
585 |
|
|
kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*& |
586 |
|
|
(-mp(i,2)-imp(i,j)/RG*dxpres(i,j))& |
587 |
|
|
+ (imp(i,j)/RG*dxpres(i,j))*(evap(i,j)/RG*dxpres(i,j))) |
588 |
|
|
else |
589 |
|
|
kappa(i,j)=((pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*& |
590 |
|
|
(-mp(i,j)-imp(i,j)/RG*dxpres(i,j))& |
591 |
|
|
+ (imp(i,j)/RG*dxpres(i,j))*(evap(i,j)/RG*dxpres(i,j))) |
592 |
|
|
endif |
593 |
|
|
else |
594 |
|
|
kappa(i,j)=1. |
595 |
|
|
endif |
596 |
|
|
ENDDO |
597 |
|
|
ENDDO |
598 |
|
|
|
599 |
|
|
DO j=klev-1,1,-1 |
600 |
|
|
DO i=1,klon |
601 |
|
|
if (abs(kappa(i,j)).lt.1.e-25) then !si denominateur nul (il peut y avoir des mp!=0) |
602 |
|
|
kappa(i,j)=1. |
603 |
|
|
if(j.eq.1) then |
604 |
|
|
qDi(i,j,it)=qDi(i,j+1,it) !orig tr(i,j,it) ! mp(1)=0 donc tout vient de la couche supérieure |
605 |
|
|
elseif(mp(i,j+1).gt.mp(i,j).and.mp(i,j+1).gt.1.e-10) then |
606 |
|
|
qDi(i,j,it)=qDi(i,j+1,it) |
607 |
|
|
elseif(mp(i,j).gt.mp(i,j+1).and.mp(i,j).gt.1.e-10) then! entrainement |
608 |
|
|
qDi(i,j,it)=(-mp(i,j+1)*(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it))/(-mp(i,j)) |
609 |
|
|
else ! si mp (i)=0 et mp(j+1)=0 |
610 |
|
|
qDi(i,j,it)=tr(i,j,it) ! orig 0. |
611 |
|
|
endif |
612 |
|
|
|
613 |
|
|
if(NO_precip(i,j)) then |
614 |
|
|
qPr(i,j,it)=0. |
615 |
|
|
else |
616 |
|
|
qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+& |
617 |
|
|
Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)& |
618 |
|
|
+imp(i,j)/RG*dxpres(i,j)*qDi(i,j,it))/& |
619 |
|
|
(pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j)) |
620 |
|
|
endif |
621 |
|
|
else ! denominateur non nul |
622 |
|
|
kappa(i,j)=1./kappa(i,j) |
623 |
|
|
! calcul de qd et qp |
624 |
|
|
!!jyg (20130119) correction pour le sommet du nuage |
625 |
|
|
!! if(j.ge.inb(i)) then !au-dessus du nuage, sommet inclu |
626 |
|
|
if(j.gt.inb(i)) then !au-dessus du nuage |
627 |
|
|
qDi(i,j,it)=tr(i,j,it) ! pas de descente => environnement = descente insaturee |
628 |
|
|
qPr(i,j,it)=0. |
629 |
|
|
|
630 |
|
|
! vvv premiere couche du modele ou mp(1)=0 ! det tout le temps vvv |
631 |
|
|
elseif(j.eq.1) then |
632 |
|
|
if(mp(i,2).gt.1.e-10) then !mp(2) non nul -> detrainement (car mp(1) = 0) !ent pas possible |
633 |
|
|
if(NO_precip(i,j)) then !pas de precip en (i) |
634 |
|
|
qDi(i,j,it)=qDi(i,j+1,it) |
635 |
|
|
qPr(i,j,it)=0. |
636 |
|
|
else |
637 |
|
|
qDi(i,j,it)=kappa(i,j)*(& |
638 |
|
|
(-evap(i,j)/RG*dxpres(i,j))*((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+& |
639 |
|
|
Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)) +& |
640 |
|
|
(pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*& |
641 |
|
|
(-mp(i,j+1)*qDi(i,j+1,it))) |
642 |
|
|
|
643 |
|
|
qPr(i,j,it)=kappa(i,j)*(& |
644 |
|
|
(-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))*& |
645 |
|
|
((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+& |
646 |
|
|
Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it))& |
647 |
|
|
+(-mp(i,j+1)*qDi(i,j+1,it)) * (imp(i,j)/RG*dxpres(i,j))) |
648 |
|
|
endif |
649 |
|
|
|
650 |
|
|
else !mp(2) nul -> plus de descente insaturee -> pluie agit sur environnement |
651 |
|
|
qDi(i,j,it)=tr(i,j,it) ! orig 0. |
652 |
|
|
if(NO_precip(i,j)) then |
653 |
|
|
qPr(i,j,it)=0. |
654 |
|
|
else |
655 |
|
|
qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+& |
656 |
|
|
Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)& |
657 |
|
|
+imp(i,j)/RG*dxpres(i,j)*tr(i,j,it))/& |
658 |
|
|
(pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j)) |
659 |
|
|
endif |
660 |
|
|
|
661 |
|
|
endif !mp(2) nul ou non |
662 |
|
|
|
663 |
|
|
! vvv (j!=1.and.j.lt.inb(i)) en-dessous du sommet nuage vvv |
664 |
|
|
else |
665 |
|
|
!------------------------------------------------------------- detrainement |
666 |
|
|
if(mp(i,j+1).gt.mp(i,j).and.mp(i,j+1).gt.1.e-10) then !mp(i,j).gt.1.e-10) then |
667 |
|
|
if(NO_precip(i,j)) then |
668 |
|
|
qDi(i,j,it)=qDi(i,j+1,it) |
669 |
|
|
qPr(i,j,it)=0. |
670 |
|
|
else |
671 |
|
|
qDi(i,j,it)=kappa(i,j)*(& |
672 |
|
|
(-evap(i,j)/RG*dxpres(i,j))*((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+& |
673 |
|
|
Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)) +& |
674 |
|
|
(pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*& |
675 |
|
|
(-mp(i,j+1)*qDi(i,j+1,it))) |
676 |
|
|
! |
677 |
|
|
qPr(i,j,it)=kappa(i,j)*(& |
678 |
|
|
(-mp(i,j+1)-imp(i,j)/RG*dxpres(i,j))*& |
679 |
|
|
((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+& |
680 |
|
|
Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it))& |
681 |
|
|
+(-mp(i,j+1)*qDi(i,j+1,it)) * (imp(i,j)/RG*dxpres(i,j))) |
682 |
|
|
endif !precip |
683 |
|
|
!------------------------------------------------------------- entrainement |
684 |
|
|
elseif(mp(i,j).gt.mp(i,j+1).and.mp(i,j).gt.1.e-10) then |
685 |
|
|
if(NO_precip(i,j)) then |
686 |
|
|
qDi(i,j,it)=(-mp(i,j+1)*(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it))/(-mp(i,j)) |
687 |
|
|
qPr(i,j,it)=0. |
688 |
|
|
else |
689 |
|
|
qDi(i,j,it)=kappa(i,j)*(& |
690 |
|
|
(-evap(i,j)/RG*dxpres(i,j))*((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+& |
691 |
|
|
Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)) +& |
692 |
|
|
(pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j))*& |
693 |
|
|
(-mp(i,j+1)*(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it))) |
694 |
|
|
! |
695 |
|
|
qPr(i,j,it)=kappa(i,j)*(& |
696 |
|
|
(-mp(i,j)-imp(i,j)/RG*dxpres(i,j))*& |
697 |
|
|
((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+& |
698 |
|
|
Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it))& |
699 |
|
|
+(-mp(i,j+1)*(qDi(i,j+1,it)-tr(i,j,it))-mp(i,j)*tr(i,j,it))*& |
700 |
|
|
(imp(i,j)/RG*dxpres(i,j))) |
701 |
|
|
endif !precip |
702 |
|
|
!------------------------------------------------------------- endif ! ent/det |
703 |
|
|
else !mp nul |
704 |
|
|
qDi(i,j,it)=tr(i,j,it) ! orig 0. |
705 |
|
|
if(NO_precip(i,j)) then |
706 |
|
|
qPr(i,j,it)=0. |
707 |
|
|
else |
708 |
|
|
qPr(i,j,it)=((pmflxr(i,j+1)+pmflxs(i,j+1))*qPr(i,j+1,it)+& |
709 |
|
|
Pa(i,j)*qPa(i,j,it)+Pm(i,j)*qMel(i,j,it)& |
710 |
|
|
+imp(i,j)/RG*dxpres(i,j)*tr(i,j,it))/& |
711 |
|
|
(pmflxr(i,j+1)+pmflxs(i,j+1)+Pa(i,j)+Pm(i,j)) |
712 |
|
|
endif |
713 |
|
|
endif ! mp nul ou non |
714 |
|
|
endif ! condition sur j |
715 |
|
|
endif ! kappa |
716 |
|
|
ENDDO |
717 |
|
|
ENDDO |
718 |
|
|
|
719 |
|
|
!! print test descente insaturee |
720 |
|
|
! DO j=klev,1,-1 |
721 |
|
|
! DO i=1,klon |
722 |
|
|
! if(it.eq.3) then |
723 |
|
|
! write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)') j,& |
724 |
|
|
!! 'zmfdam',zmfdam(i,j,it),'zmfpsi',zmfphi2(i,j,it),& |
725 |
|
|
! 'zmfdam+zmfpsi',zmfdam(i,j,it)+zmfphi2(i,j,it),'qpmMint',qpmMint(i,j,it),& |
726 |
|
|
! 'Pm',Pm(i,j),'Mint',Mint(i,j),& |
727 |
|
|
!! 'zmfa',zmfa(i,j,it),'zmfp',zmfp(i,j,it),& |
728 |
|
|
! 'zmfdam',zmfdam(i,j,it),'zmfpsi',zmfphi2(i,j,it),'zmfd1a',zmfd1a(i,j,it) |
729 |
|
|
!! 'Pa',Pa(i,j),'eplaMm',eplaMm(i,j) |
730 |
|
|
!! 'zmfd1a=da1*qa',zmfd1a(i,j,it),'Pa*qPa',wdtrainA(i,j)*qPa(i,j,it),'da1',d1a(i,j) |
731 |
|
|
! endif |
732 |
|
|
! ENDDO |
733 |
|
|
! ENDDO |
734 |
|
|
|
735 |
|
|
|
736 |
|
|
! =================================================== |
737 |
|
|
! calcul final des tendances |
738 |
|
|
! =================================================== |
739 |
|
|
|
740 |
|
|
DO k=klev-1,1,-1 |
741 |
|
|
DO i=1, klon |
742 |
|
|
! transport |
743 |
|
|
tdcvMA=zmfd(i,k,it)+zmfu(i,k,it)+zmfa(i,k,it)+zmfp(i,k,it) ! double comptage des downdraft insatures |
744 |
|
|
trsptrac=zmfu(i,k,it)+zmfa(i,k,it)+zmfp(i,k,it) |
745 |
|
|
! lessivage courants satures |
746 |
|
|
!JE scavtrac=-ccntrAA_coef*zmfd1a(i,k,it)& |
747 |
|
|
!JE -zmfphi2(i,k,it)*ccntrENV_coef& |
748 |
|
|
!JE -zmfdam(i,k,it)*ccntrAA_coef |
749 |
|
|
scavtrac=-ccntrAA_coef3d(i,k)*zmfd1a(i,k,it)& |
750 |
|
|
-zmfphi2(i,k,it)*ccntrENV_coef3d(i,k)& |
751 |
|
|
-zmfdam(i,k,it)*ccntrAA_coef3d(i,k) |
752 |
|
|
! lessivage courants insatures |
753 |
|
|
if(k.le.inb(i).and.k.gt.1) then ! tendances dans le nuage |
754 |
|
|
!------------------------------------------------------------- detrainement |
755 |
|
|
if(mp(i,k+1).gt.mp(i,k).and.mp(i,k+1).gt.1.e-10) then |
756 |
|
|
uscavtrac= (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))& |
757 |
|
|
+ mp(i,k)*(tr(i,k-1,it)-tr(i,k,it)) |
758 |
|
|
! |
759 |
|
|
! if(it.eq.3) write(*,'(I2,1X,a,5X,e20.12,82X,a,e20.12)')k,' det incloud',& |
760 |
|
|
! (-mp(i,k)+mp(i,k+1))*(qDi(i,k,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k)+& |
761 |
|
|
! mp(i,k)*(tr(i,k-1,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k),& |
762 |
|
|
! 'mp',mp(i,k) |
763 |
|
|
!------------------------------------------------------------- entrainement |
764 |
|
|
elseif(mp(i,k).gt.mp(i,k+1).and.mp(i,k).gt.1.e-10) then |
765 |
|
|
uscavtrac= mp(i,k)*(tr(i,k-1,it)-tr(i,k,it)) |
766 |
|
|
! |
767 |
|
|
! if(it.eq.3) write(*,'(I2,1X,a,5X,e20.12,82X,a,e20.12)')k,' ent incloud',uscavtrac*pdtimeRG/dxpres(i,k), 'mp',mp(i,k) |
768 |
|
|
!=!------------------------------------------------------------- end ent/det |
769 |
|
|
else ! mp(i,k+1)=0. et mp(i,k)=0. pluie directement sur l environnement |
770 |
|
|
|
771 |
|
|
if(NO_precip(i,k)) then |
772 |
|
|
uscavtrac=0. |
773 |
|
|
! if(it.eq.3) write(*,'(I2,1X,a,e20.12,82X,a,e20.12)')k,' no P ent incloud',uscavtrac*pdtimeRG/dxpres(i,k), 'mp',mp(i,k) |
774 |
|
|
else |
775 |
|
|
uscavtrac=-imp(i,k)*tr(i,k,it)*dxpres(i,k)/RG+evap(i,k)*qPr(i,k,it)*dxpres(i,k)/RG |
776 |
|
|
! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac*pdtimeRG/dxpres(i,k), 'mp',mp(i,k) |
777 |
|
|
!!JE adds |
778 |
|
|
! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'imp',imp(i,k) |
779 |
|
|
! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'tr',tr(i,k,it) |
780 |
|
|
! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'evap',evap(i,k) |
781 |
|
|
! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'qPr',qPr(i,k,it) |
782 |
|
|
! if(it.eq.3) write(*,'(I2,1X,a,3X,e20.12,82X,a,e20.12)')k,' P env incloud',uscavtrac, 'dxpres',dxpres(i,k) |
783 |
|
|
!!Je end |
784 |
|
|
|
785 |
|
|
endif |
786 |
|
|
endif ! mp/det/ent |
787 |
|
|
!------------------------------------------------------------- premiere couche |
788 |
|
|
elseif(k.eq.1) then ! mp(1)=0. |
789 |
|
|
if(mp(i,2).gt.1.e-10) then !detrainement |
790 |
|
|
uscavtrac= (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it)) !& |
791 |
|
|
! + mp(i,2)*(0.-tr(i,k,it)) |
792 |
|
|
! |
793 |
|
|
! if(it.eq.3) write(*,'(I2,1X,a,e20.12,84X,a,e20.12)')k,' 1 det',& |
794 |
|
|
! (-0.+mp(i,2))*(qDi(i,k,it)-tr(i,k,it))*pdtimeRG/dxpres(i,k)+& |
795 |
|
|
! mp(i,2)*(0.-tr(i,k,it))*pdtimeRG/dxpres(i,k),& |
796 |
|
|
! 'mp',mp(i,k) |
797 |
|
|
else ! mp(2) = 0 = mp(1) pas de descente insaturee, rien ne se passe s'il ne pleut pas, sinon pluie->env |
798 |
|
|
if(NO_precip(i,1)) then |
799 |
|
|
uscavtrac=0. |
800 |
|
|
else |
801 |
|
|
uscavtrac=-imp(i,k)*tr(i,k,it)*dxpres(i,k)/RG+evap(i,k)*qPr(i,k,it)*dxpres(i,k)/RG |
802 |
|
|
endif |
803 |
|
|
! if(it.eq.3) write(*,'(I2,1X,a,2X,e20.12,82X,a,e20.12)')k,'1 P env incloud',uscavtrac*pdtimeRG/dxpres(i,k), 'mp',mp(i,k) |
804 |
|
|
endif |
805 |
|
|
|
806 |
|
|
else ! k > INB au-dessus du nuage |
807 |
|
|
uscavtrac=0. |
808 |
|
|
endif |
809 |
|
|
|
810 |
|
|
! ===== tendances finales ====== |
811 |
|
|
trsptd(i,k,it)=trsptrac*pdtimeRG/dxpres(i,k) ! td transport sans eau dans courants satures |
812 |
|
|
dtrSscav(i,k,it)=scavtrac*pdtimeRG/dxpres(i,k) ! td du lessivage dans courants satures |
813 |
|
|
dtrUscav(i,k,it)=uscavtrac*pdtimeRG/dxpres(i,k) ! td courant insat |
814 |
|
|
dtrsat(i,k,it)=(trsptrac+scavtrac)*pdtimeRG/dxpres(i,k) ! td courant sat |
815 |
|
|
dtrcv(i,k,it)=(trsptrac+scavtrac+uscavtrac)*pdtimeRG/dxpres(i,k)!dtrsat(i,k,it)+dtrUscav(i,k,it) td conv |
816 |
|
|
!!!!!! |
817 |
|
|
dtrcvMA(i,k,it)=tdcvMA*pdtimeRG/dxpres(i,k) ! MA tendance convection |
818 |
|
|
ENDDO |
819 |
|
|
ENDDO |
820 |
|
|
|
821 |
|
|
! test de conservation du traceur |
822 |
|
|
!print*,"_____________________________________________________________" |
823 |
|
|
!print*," " |
824 |
|
|
! conserv=0. |
825 |
|
|
! conservMA=0. |
826 |
|
|
! DO k= klev-1,1,-1 |
827 |
|
|
! DO i=1, klon |
828 |
|
|
! conserv=conserv+dtrcv(i,k,it)* & |
829 |
|
|
! (paprs(i,k)-paprs(i,k+1))/RG |
830 |
|
|
! conservMA=conservMA+dtrcvMA(i,k,it)* & |
831 |
|
|
! (paprs(i,k)-paprs(i,k+1))/RG |
832 |
|
|
! |
833 |
|
|
! if(it.eq.3) write(*,'(I2,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12,2X,a,e20.12)') k,& |
834 |
|
|
! 'MA td ',dtrcvMA(i,k,it)*dxpres(i,k)/RG,& |
835 |
|
|
! ' td',dtrcv(i,k,it)*dxpres(i,k)/RG,' conservMA ',conservMA,'conserv ',conserv |
836 |
|
|
!! |
837 |
|
|
! ENDDO |
838 |
|
|
! ENDDO |
839 |
|
|
! if(it.eq.3) print *,'it',it,'conserv ',conserv,'conservMA ',conservMA |
840 |
|
|
|
841 |
|
|
END SUBROUTINE cvltr_spl |