Context Navigation

phytrac_mod.F90 @ 4250

Last change on this file since 4250 was 4170, checked in by dcugnet, 2 years ago
The variable "types_trac" is the equivalent of "type_trac" in case multiple sections must be read and used in "tracer.def" file. Tests on the "type_trac" were replaced with tests on the vector "types_trac". Most of the time, there are two components: 'lmdz' and a second one. The later has priority on 'lmdz' and must be used for the tests. For more components, care must be taken to execute specific parts of the code on the right tracers ; the tracers(:)%component has been created in that respect.
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:keywords set to `Id`
File size: 45.5 KB

Line
1	!$Id: phytrac_mod.F90 4170 2022-06-16 18:16:59Z oboucher $
2	MODULE phytrac_mod
3	!=================================================================================
4	! Interface between the LMDZ physical package and tracer computation.
5	! Chemistry modules (INCA, Reprobus or the more specific traclmdz routine)
6	! are called from phytrac.
7	!
8	!======================================================================
9	! Auteur(s) FH
10	! Objet: Moniteur general des tendances traceurs
11	!
12	! iflag_vdf_trac : Options for activating transport by vertical diffusion :
13	! 1. notmal
14	! 0. emission is injected in the first layer only, without diffusion
15	! -1 no emission & no diffusion
16	! Modification 2013/07/22 : transformed into a module to pass tendencies to
17	! physics outputs. Additional keys for controling activation of sub processes.
18	! Modification R. Pilon 10 octobre 2012 large scale scavenging incloud_scav + bc_scav
19	! Modification R. Pilon 01 janvier 2012 transport+scavenging KE scheme : cvltr
20	!=================================================================================
21
22	!
23	! Tracer tendencies, for outputs
24	!-------------------------------
25	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_cl ! Td couche limite/traceur
26	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_dec !RomP
27	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_cv ! Td convection/traceur
28	! RomP >>>
29	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_insc
30	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_bcscav
31	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_evapls
32	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_ls
33	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_trsp
34	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_sscav
35	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_sat
36	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_uscav
37	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: qPr,qDi ! concentration tra dans pluie,air descente insaturee
38	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: qPa,qMel
39	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: qTrdi,dtrcvMA ! conc traceur descente air insaturee et td convective MA
40	! RomP <<<
41	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_th ! Td thermique
42	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_lessi_impa ! Td du lessivage par impaction
43	REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: d_tr_lessi_nucl ! Td du lessivage par nucleation
44	REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: qPrls !jyg: concentration tra dans pluie LS a la surf.
45	REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: d_tr_dry ! Td depot sec/traceur (1st layer),ALLOCATABLE,SAVE jyg
46	REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: flux_tr_dry ! depot sec/traceur (surface),ALLOCATABLE,SAVE jyg
47
48	!$OMP THREADPRIVATE(qPa,qMel,qTrdi,dtrcvMA,d_tr_th,d_tr_lessi_impa,d_tr_lessi_nucl)
49	!$OMP THREADPRIVATE(d_tr_trsp,d_tr_sscav,d_tr_sat,d_tr_uscav,qPr,qDi)
50	!$OMP THREADPRIVATE(d_tr_insc,d_tr_bcscav,d_tr_evapls,d_tr_ls,qPrls)
51	!$OMP THREADPRIVATE(d_tr_cl,d_tr_dry,flux_tr_dry,d_tr_dec,d_tr_cv)
52
53
54	CONTAINS
55
56	SUBROUTINE phytrac_init()
57	USE dimphy
58	USE infotrac_phy, ONLY: nbtr, types_trac
59	USE tracco2i_mod, ONLY: tracco2i_init
60	IMPLICIT NONE
61
62	ALLOCATE(d_tr_cl(klon,klev,nbtr),d_tr_dry(klon,nbtr))
63	ALLOCATE(flux_tr_dry(klon,nbtr),d_tr_dec(klon,klev,nbtr),d_tr_cv(klon,klev,nbtr))
64	ALLOCATE(d_tr_insc(klon,klev,nbtr),d_tr_bcscav(klon,klev,nbtr))
65	ALLOCATE(d_tr_evapls(klon,klev,nbtr),d_tr_ls(klon,klev,nbtr))
66	ALLOCATE(qPrls(klon,nbtr),d_tr_trsp(klon,klev,nbtr))
67	ALLOCATE(d_tr_sscav(klon,klev,nbtr),d_tr_sat(klon,klev,nbtr))
68	ALLOCATE(d_tr_uscav(klon,klev,nbtr),qPr(klon,klev,nbtr),qDi(klon,klev,nbtr))
69	ALLOCATE(qPa(klon,klev,nbtr),qMel(klon,klev,nbtr))
70	ALLOCATE(qTrdi(klon,klev,nbtr),dtrcvMA(klon,klev,nbtr))
71	ALLOCATE(d_tr_th(klon,klev,nbtr))
72	ALLOCATE(d_tr_lessi_impa(klon,klev,nbtr),d_tr_lessi_nucl(klon,klev,nbtr))
73
74
75
76	!===============================================================================
77	! -- Do specific treatment according to chemestry model or local LMDZ tracers
78	!
79	!===============================================================================
80	! -- CO2 interactif --
81	IF(ANY(types_trac == 'co2i') .OR. ANY(types_trac == 'inco')) CALL tracco2i_init()
82
83
84	END SUBROUTINE phytrac_init
85
86	SUBROUTINE phytrac( &
87	nstep, julien, gmtime, debutphy, &
88	lafin, pdtphys, u, v, t_seri, &
89	paprs, pplay, pmfu, pmfd, &
90	pen_u, pde_u, pen_d, pde_d, &
91	cdragh, coefh, fm_therm, entr_therm, &
92	yu1, yv1, ftsol, pctsrf, &
93	ustar, u10m, v10m, &
94	wstar, ale_bl, ale_wake, &
95	xlat, xlon, &
96	frac_impa,frac_nucl,beta_fisrt,beta_v1, &
97	presnivs, pphis, pphi, albsol, &
98	sh, ch, rh, cldfra, rneb, &
99	diafra, cldliq, itop_con, ibas_con, &
100	pmflxr, pmflxs, prfl, psfl, &
101	da, phi, mp, upwd, &
102	phi2, d1a, dam, sij, wght_cvfd, & ! RomP +RL
103	wdtrainA, wdtrainM, sigd, clw, elij, & ! RomP
104	evap, ep, epmlmMm, eplaMm, & ! RomP
105	dnwd, aerosol_couple, flxmass_w, &
106	tau_aero, piz_aero, cg_aero, ccm, &
107	rfname, &
108	d_tr_dyn, & ! RomP
109	tr_seri, init_source)
110	!
111	!======================================================================
112	! Auteur(s) FH
113	! Objet: Moniteur general des tendances traceurs
114	! Modification R. Pilon 01 janvier 2012 transport+scavenging KE scheme : cvltr
115	! Modification R. Pilon 10 octobre 2012 large scale scavenging incloud_scav + bc_scav
116	!======================================================================
117
118	USE ioipsl
119	USE phys_cal_mod, only : hour
120	USE dimphy
121	USE infotrac_phy, ONLY: nbtr, nqCO2, types_trac, type_trac, conv_flg, pbl_flg
122	USE strings_mod, ONLY: int2str
123	USE mod_grid_phy_lmdz
124	USE mod_phys_lmdz_para
125	USE iophy
126	USE traclmdz_mod
127	USE tracinca_mod
128	USE tracreprobus_mod
129	USE indice_sol_mod
130	USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root
131	USE print_control_mod, ONLY: lunout
132	USE aero_mod, ONLY : naero_grp
133
134	USE tracco2i_mod
135
136	#ifdef CPP_StratAer
137	USE traccoag_mod
138	USE phys_local_var_mod, ONLY: mdw
139	USE phys_local_var_mod, ONLY: budg_dep_dry_ocs, budg_dep_wet_ocs
140	USE phys_local_var_mod, ONLY: budg_dep_dry_so2, budg_dep_wet_so2
141	USE phys_local_var_mod, ONLY: budg_dep_dry_h2so4, budg_dep_wet_h2so4
142	USE phys_local_var_mod, ONLY: budg_dep_dry_part, budg_dep_wet_part
143	USE infotrac_phy, ONLY: nbtr_sulgas, id_OCS_strat, id_SO2_strat, id_H2SO4_strat
144	USE aerophys
145	#endif
146
147	IMPLICIT NONE
148
149	INCLUDE "YOMCST.h"
150	INCLUDE "clesphys.h"
151	!==========================================================================
152	! -- ARGUMENT DESCRIPTION --
153	!==========================================================================
154
155	! Input arguments
156	!----------------
157	!Configuration grille,temps:
158	INTEGER,INTENT(IN) :: nstep ! Appel physique
159	INTEGER,INTENT(IN) :: julien ! Jour julien
160	REAL,INTENT(IN) :: gmtime ! Heure courante
161	REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde)
162	LOGICAL,INTENT(IN) :: debutphy ! le flag de l'initialisation de la physique
163	LOGICAL,INTENT(IN) :: lafin ! le flag de la fin de la physique
164
165	REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point
166	REAL,DIMENSION(klon),INTENT(IN) :: xlon ! longitudes pour chaque point
167	!
168	!Physique:
169	!--------
170	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
171	REAL,DIMENSION(klon,klev),INTENT(IN) :: u ! variable not used
172	REAL,DIMENSION(klon,klev),INTENT(IN) :: v ! variable not used
173	REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique
174	REAL,DIMENSION(klon,klev),INTENT(IN) :: rh ! humidite relative
175	REAL,DIMENSION(klon,klev),INTENT(IN) :: ch ! eau liquide (+ glace si le traceur existe)
176	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa)
177	REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa)
178	REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi ! geopotentiel
179	REAL,DIMENSION(klon),INTENT(IN) :: pphis
180	REAL,DIMENSION(klev),INTENT(IN) :: presnivs
181	REAL,DIMENSION(klon,klev),INTENT(IN) :: cldliq ! eau liquide nuageuse
182	REAL,DIMENSION(klon,klev),INTENT(IN) :: cldfra ! fraction nuageuse (tous les nuages)
183	REAL,DIMENSION(klon,klev),INTENT(IN) :: diafra ! fraction nuageuse (convection ou stratus artificiels)
184	REAL,DIMENSION(klon,klev),INTENT(IN) :: rneb ! fraction nuageuse (grande echelle)
185	!
186	REAL :: ql_incl ! contenu en eau liquide nuageuse dans le nuage ! ql_incl=oliq/rneb
187	REAL,DIMENSION(klon,klev),INTENT(IN) :: beta_fisrt ! taux de conversion de l'eau cond (de fisrtilp)
188	REAL,DIMENSION(klon,klev),INTENT(out) :: beta_v1 ! -- (originale version)
189
190	!
191	INTEGER,DIMENSION(klon),INTENT(IN) :: itop_con
192	INTEGER,DIMENSION(klon),INTENT(IN) :: ibas_con
193	REAL,DIMENSION(klon),INTENT(IN) :: albsol ! albedo surface
194	!
195	!Dynamique
196	!--------
197	REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: d_tr_dyn
198	!
199	!Convection:
200	!----------
201	REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfu ! flux de masse dans le panache montant
202	REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfd ! flux de masse dans le panache descendant
203	REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_u ! flux entraine dans le panache montant
204	REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_u ! flux detraine dans le panache montant
205	REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_d ! flux entraine dans le panache descendant
206	REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_d ! flux detraine dans le panache descendant
207
208	!...Tiedke
209	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: pmflxr, pmflxs ! Flux precipitant de pluie, neige aux interfaces [convection]
210	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: prfl, psfl ! Flux precipitant de pluie, neige aux interfaces [large-scale]
211
212	LOGICAL,INTENT(IN) :: aerosol_couple
213	REAL,DIMENSION(klon,klev),INTENT(IN) :: flxmass_w
214	REAL,DIMENSION(klon,klev,naero_grp,2),INTENT(IN) :: tau_aero
215	REAL,DIMENSION(klon,klev,naero_grp,2),INTENT(IN) :: piz_aero
216	REAL,DIMENSION(klon,klev,naero_grp,2),INTENT(IN) :: cg_aero
217	CHARACTER(len=4),DIMENSION(naero_grp),INTENT(IN) :: rfname
218	REAL,DIMENSION(klon,klev,2),INTENT(IN) :: ccm
219	!... K.Emanuel
220	REAL,DIMENSION(klon,klev),INTENT(IN) :: da
221	REAL,DIMENSION(klon,klev,klev),INTENT(IN):: phi
222	! RomP >>>
223	REAL,DIMENSION(klon,klev),INTENT(IN) :: d1a,dam
224	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi2
225	!
226	REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtrainA
227	REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtrainM
228	REAL,DIMENSION(klon),INTENT(IN) :: sigd
229	! ---- RomP flux entraine, detraine et precipitant kerry Emanuel
230	REAL,DIMENSION(klon,klev),INTENT(IN) :: evap
231	REAL,DIMENSION(klon,klev),INTENT(IN) :: ep
232	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: sij
233	REAL,DIMENSION(klon,klev),INTENT(IN) :: wght_cvfd !RL
234	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: elij
235	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: epmlmMm
236	REAL,DIMENSION(klon,klev),INTENT(IN) :: eplaMm
237	REAL,DIMENSION(klon,klev),INTENT(IN) :: clw
238	! RomP <<<
239
240	!
241	REAL,DIMENSION(klon,klev),INTENT(IN) :: mp
242	REAL,DIMENSION(klon,klev),INTENT(IN) :: upwd ! saturated updraft mass flux
243	REAL,DIMENSION(klon,klev),INTENT(IN) :: dnwd ! saturated downdraft mass flux
244	!
245	!Thermiques:
246	!----------
247	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: fm_therm
248	REAL,DIMENSION(klon,klev),INTENT(IN) :: entr_therm
249	!
250	!Couche limite:
251	!--------------
252	!
253	REAL,DIMENSION(:),INTENT(IN) :: cdragh ! (klon) coeff drag pour T et Q
254	REAL,DIMENSION(:,:),INTENT(IN) :: coefh ! (klon,klev) coeff melange CL (m**2/s)
255	REAL,DIMENSION(:),INTENT(IN) :: ustar,u10m,v10m ! (klon) u* & vent a 10m (m/s)
256	REAL,DIMENSION(:),INTENT(IN) :: wstar,ale_bl,ale_wake ! (klon) w* and Avail. Lifting Ener.
257	REAL,DIMENSION(:),INTENT(IN) :: yu1 ! (klon) vents au premier niveau
258	REAL,DIMENSION(:),INTENT(IN) :: yv1 ! (klon) vents au premier niveau
259	!
260	!Lessivage:
261	!----------
262	!
263	REAL, DIMENSION(:), ALLOCATABLE, SAVE :: ccntrAA
264	REAL, DIMENSION(:), ALLOCATABLE, SAVE :: ccntrENV
265	REAL, DIMENSION(:), ALLOCATABLE, SAVE :: coefcoli
266	LOGICAL, DIMENSION(:), ALLOCATABLE, SAVE :: flag_cvltr
267	!$OMP THREADPRIVATE(ccntrAA,ccntrENV,coefcoli,flag_cvltr)
268	REAL, DIMENSION(klon,klev) :: ccntrAA_3d
269	REAL, DIMENSION(klon,klev) :: ccntrENV_3d
270	REAL, DIMENSION(klon,klev) :: coefcoli_3d
271	!
272	! pour le ON-LINE
273	!
274	REAL,DIMENSION(klon,klev),INTENT(IN) :: frac_impa ! fraction d'aerosols non impactes
275	REAL,DIMENSION(klon,klev),INTENT(IN) :: frac_nucl ! fraction d'aerosols non nuclees
276
277	! Arguments necessaires pour les sources et puits de traceur:
278	REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol ! Temperature du sol (surf)(Kelvin)
279	REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol (nature du sol)
280
281	#ifdef CPP_StratAer
282	REAL,DIMENSION(klon) :: v_dep_dry !dry deposition velocity of stratospheric sulfate in m/s
283	#endif
284	! Output argument
285	!----------------
286	REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
287	REAL,DIMENSION(klon,klev) :: sourceBE
288	REAL,DIMENSION(klon,nbtr), INTENT(IN) :: init_source
289
290	!=======================================================================================
291	! -- LOCAL VARIABLES --
292	!=======================================================================================
293
294	INTEGER :: i, k, it
295	INTEGER :: nsplit
296
297	!Sources et Reservoirs de traceurs (ex:Radon):
298	!--------------------------------------------
299	!
300	REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: source ! a voir lorsque le flux de surface est prescrit
301	!$OMP THREADPRIVATE(source)
302
303	!
304	!Entrees/Sorties:
305	!---------------
306	INTEGER :: iiq, ierr
307	INTEGER :: nhori, nvert
308	REAL :: zsto, zout, zjulian
309	INTEGER,SAVE :: nid_tra ! pointe vers le fichier histrac.nc
310	!$OMP THREADPRIVATE(nid_tra)
311	REAL,DIMENSION(klon) :: zx_tmp_fi2d ! variable temporaire grille physique
312	INTEGER :: itau_w ! pas de temps ecriture = nstep + itau_phy
313	LOGICAL,PARAMETER :: ok_sync=.TRUE.
314	!
315	! Nature du traceur
316	!------------------
317	LOGICAL,DIMENSION(:),ALLOCATABLE,SAVE :: aerosol ! aerosol(it) = true => aerosol => lessivage
318	!$OMP THREADPRIVATE(aerosol) ! aerosol(it) = false => gaz
319	REAL,DIMENSION(klon,klev) :: delp ! epaisseur de couche (Pa)
320	!
321	! Tendances de traceurs (Td) et flux de traceurs:
322	!------------------------
323	REAL,DIMENSION(klon,klev) :: d_tr ! Td dans l'atmosphere
324	REAL,DIMENSION(klon,klev) :: Mint
325	REAL,DIMENSION(klon,klev,nbtr) :: zmfd1a
326	REAL,DIMENSION(klon,klev,nbtr) :: zmfdam
327	REAL,DIMENSION(klon,klev,nbtr) :: zmfphi2
328	! Physique
329	!----------
330	REAL,DIMENSION(klon,klev,nbtr) :: flestottr ! flux de lessivage dans chaque couche
331	REAL,DIMENSION(klon,klev) :: zmasse ! densite atmospherique Kg/m2
332	REAL,DIMENSION(klon,klev) :: ztra_th
333	!PhH
334	REAL,DIMENSION(klon,klev) :: zrho
335	REAL,DIMENSION(klon,klev) :: zdz
336	REAL :: evaplsc,dx,beta ! variable pour lessivage Genthon
337	REAL,DIMENSION(klon) :: his_dh ! ---
338	! in-cloud scav variables
339	REAL :: ql_incloud_ref ! ref value of in-cloud condensed water content
340
341	!Controles:
342	!---------
343	INTEGER,SAVE :: iflag_vdf_trac,iflag_con_trac,iflag_the_trac
344	INTEGER,SAVE :: iflag_con_trac_omp, iflag_vdf_trac_omp,iflag_the_trac_omp
345	!$OMP THREADPRIVATE(iflag_vdf_trac,iflag_con_trac,iflag_the_trac)
346
347	LOGICAL,SAVE :: lessivage
348	!$OMP THREADPRIVATE(lessivage)
349
350	!RomP >>>
351	INTEGER,SAVE :: iflag_lscav_omp,iflag_lscav
352	REAL, SAVE :: ccntrAA_in,ccntrAA_omp
353	REAL, SAVE :: ccntrENV_in,ccntrENV_omp
354	REAL, SAVE :: coefcoli_in,coefcoli_omp
355
356	LOGICAL,SAVE :: convscav_omp,convscav
357	!$OMP THREADPRIVATE(iflag_lscav)
358	!$OMP THREADPRIVATE(ccntrAA_in,ccntrENV_in,coefcoli_in)
359	!$OMP THREADPRIVATE(convscav)
360	!RomP <<<
361	!######################################################################
362	! -- INITIALIZATION --
363	!######################################################################
364
365	DO k=1,klev
366	DO i=1,klon
367	sourceBE(i,k)=0.
368	Mint(i,k)=0.
369	zrho(i,k)=0.
370	zdz(i,k)=0.
371	END DO
372	END DO
373
374	DO it=1, nbtr
375	DO k=1,klev
376	DO i=1,klon
377	d_tr_insc(i,k,it)=0.
378	d_tr_bcscav(i,k,it)=0.
379	d_tr_evapls(i,k,it)=0.
380	d_tr_ls(i,k,it)=0.
381	d_tr_cv(i,k,it)=0.
382	d_tr_cl(i,k,it)=0.
383	d_tr_trsp(i,k,it)=0.
384	d_tr_sscav(i,k,it)=0.
385	d_tr_sat(i,k,it)=0.
386	d_tr_uscav(i,k,it)=0.
387	d_tr_lessi_impa(i,k,it)=0.
388	d_tr_lessi_nucl(i,k,it)=0.
389	qDi(i,k,it)=0.
390	qPr(i,k,it)=0.
391	qPa(i,k,it)=0.
392	qMel(i,k,it)=0.
393	qTrdi(i,k,it)=0.
394	dtrcvMA(i,k,it)=0.
395	zmfd1a(i,k,it)=0.
396	zmfdam(i,k,it)=0.
397	zmfphi2(i,k,it)=0.
398	END DO
399	END DO
400	END DO
401
402	DO it=1, nbtr
403	DO i=1,klon
404	d_tr_dry(i,it)=0.
405	flux_tr_dry(i,it)=0.
406	END DO
407	END DO
408
409	DO k = 1, klev
410	DO i = 1, klon
411	delp(i,k) = paprs(i,k)-paprs(i,k+1)
412	END DO
413	END DO
414
415	IF (debutphy) THEN
416	!!jyg
417	!$OMP BARRIER
418	ecrit_tra=86400. ! frequence de stokage en dur
419	! obsolete car remplace par des ecritures dans phys_output_write
420	!RomP >>>
421	!
422	!Config Key = convscav
423	!Config Desc = Convective scavenging switch: 0=off, 1=on.
424	!Config Def = .FALSE.
425	!Config Help =
426	!
427	!$OMP MASTER
428	convscav_omp=.FALSE.
429	call getin('convscav', convscav_omp)
430	iflag_vdf_trac_omp=1
431	call getin('iflag_vdf_trac', iflag_vdf_trac_omp)
432	iflag_con_trac_omp=1
433	call getin('iflag_con_trac', iflag_con_trac_omp)
434	iflag_the_trac_omp=1
435	call getin('iflag_the_trac', iflag_the_trac_omp)
436	!$OMP END MASTER
437	!$OMP BARRIER
438	convscav=convscav_omp
439	iflag_vdf_trac=iflag_vdf_trac_omp
440	iflag_con_trac=iflag_con_trac_omp
441	iflag_the_trac=iflag_the_trac_omp
442	write(lunout,*) 'phytrac passage dans routine conv avec lessivage', convscav
443	!
444	!Config Key = iflag_lscav
445	!Config Desc = Large scale scavenging parametrization: 0=none, 1=old(Genthon92),
446	! 2=1+PHeinrich, 3=Reddy_Boucher2004, 4=3+RPilon.
447	!Config Def = 1
448	!Config Help =
449	!
450	!$OMP MASTER
451	iflag_lscav_omp=1
452	call getin('iflag_lscav', iflag_lscav_omp)
453	ccntrAA_omp=1
454	ccntrENV_omp=1.
455	coefcoli_omp=0.001
456	call getin('ccntrAA', ccntrAA_omp)
457	call getin('ccntrENV', ccntrENV_omp)
458	call getin('coefcoli', coefcoli_omp)
459	!$OMP END MASTER
460	!$OMP BARRIER
461	iflag_lscav=iflag_lscav_omp
462	ccntrAA_in=ccntrAA_omp
463	ccntrENV_in=ccntrENV_omp
464	coefcoli_in=coefcoli_omp
465	!
466	SELECT CASE(iflag_lscav)
467	CASE(0)
468	WRITE(lunout,*) 'Large scale scavenging: none'
469	CASE(1)
470	WRITE(lunout,*) 'Large scale scavenging: C. Genthon, Tellus(1992), 44B, 371-389'
471	CASE(2)
472	WRITE(lunout,*) 'Large scale scavenging: C. Genthon, modified P. Heinrich'
473	CASE(3)
474	WRITE(lunout,*) 'Large scale scavenging: M. Shekkar Reddy and O. Boucher, JGR(2004), 109, D14202'
475	CASE(4)
476	WRITE(lunout,*) 'Large scale scavenging: Reddy and Boucher, modified R. Pilon'
477	END SELECT
478	!RomP <<<
479	WRITE(,) 'FIRST TIME IN PHYTRAC : pdtphys(sec) = ',pdtphys,'ecrit_tra (sec) = ',ecrit_tra
480	ALLOCATE( source(klon,nbtr), stat=ierr)
481	IF (ierr /= 0) CALL abort_physic('phytrac', 'pb in allocation 1',1)
482
483	ALLOCATE( aerosol(nbtr), stat=ierr)
484	IF (ierr /= 0) CALL abort_physic('phytrac', 'pb in allocation 2',1)
485
486
487	! Initialize module for specific tracers
488	IF(ANY(types_trac == 'inca')) THEN
489	source(:,:)=init_source(:,:)
490	CALL tracinca_init(aerosol,lessivage)
491	ELSE IF(ANY(types_trac == 'repr')) THEN
492	source(:,:)=0.
493	ELSE IF(ANY(types_trac == 'co2i')) THEN
494	source(:,:)=0.
495	lessivage = .FALSE.
496	aerosol(:) = .FALSE.
497	pbl_flg(:) = 1
498	iflag_the_trac= 1
499	iflag_vdf_trac= 1
500	iflag_con_trac= 1
501	ELSE IF(ANY(types_trac == 'inco')) THEN
502	source(:,1:nqCO2) = 0. ! from CO2i ModThL
503	source(:,nqCO2+1:nbtr)=init_source(:,:) ! from INCA ModThL
504	aerosol(1:nqCO2) = .FALSE. ! from CO2i ModThL
505	CALL tracinca_init(aerosol(nqCO2+1:nbtr),lessivage) ! from INCA ModThL
506	pbl_flg(1:nqCO2) = 1 ! From CO2i ModThL
507	iflag_the_trac = 1 ! From CO2i
508	iflag_vdf_trac = 1 ! From CO2i
509	iflag_con_trac = 1 ! From CO2i
510	#ifdef CPP_StratAer
511	ELSE IF(ANY(types_trac == 'coag')) THEN
512	source(:,:)=0.
513	DO it= 1, nbtr_sulgas
514	aerosol(it)=.FALSE.
515	IF (it==id_H2SO4_strat) aerosol(it)=.TRUE.
516	ENDDO
517	DO it= nbtr_sulgas+1, nbtr
518	aerosol(it)=.TRUE.
519	ENDDO
520	#endif
521	ELSE IF(ANY(types_trac == 'lmdz')) THEN
522	CALL traclmdz_init(pctsrf,xlat,xlon,ftsol,tr_seri,t_seri,pplay,sh,pdtphys,aerosol,lessivage)
523	END IF
524
525	!
526	!--initialising coefficients for scavenging in the case of NP
527	!
528	ALLOCATE(flag_cvltr(nbtr))
529	IF (iflag_con.EQ.3) THEN
530	!
531	ALLOCATE(ccntrAA(nbtr))
532	ALLOCATE(ccntrENV(nbtr))
533	ALLOCATE(coefcoli(nbtr))
534	!
535	DO it=1, nbtr
536	IF(ANY(types_trac == 'repr')) THEN
537	flag_cvltr(it)=.FALSE.
538	ELSE IF(ANY(types_trac == 'inca')) THEN
539	! IF ((it.EQ.id_Rn222) .OR. ((it.GE.id_SO2) .AND. (it.LE.id_NH3)) ) THEN
540	! !--gas-phase species
541	! flag_cvltr(it)=.FALSE.
542	!
543	! ELSEIF ( (it.GE.id_CIDUSTM) .AND. (it.LE.id_AIN) ) THEN
544	! !--insoluble aerosol species
545	! flag_cvltr(it)=.TRUE.
546	! ccntrAA(it)=0.7
547	! ccntrENV(it)=0.7
548	! coefcoli(it)=0.001
549	! ELSEIF ( (it.EQ.id_Pb210) .OR. ((it.GE.id_CSSSM) .AND. (it.LE.id_SSN))) THEN
550	! !--soluble aerosol species
551	! flag_cvltr(it)=.TRUE.
552	! ccntrAA(it)=0.9
553	! ccntrENV(it)=0.9
554	! coefcoli(it)=0.001
555	! ELSE
556	! WRITE(lunout,*) 'pb it=', it
557	! CALL abort_physic('phytrac','pb it scavenging',1)
558	! ENDIF
559	!--test OB
560	!--for now we do not scavenge in cvltr
561	flag_cvltr(it)=.FALSE.
562	ELSE IF(ANY(types_trac == 'co2i')) THEN
563	!--co2 tracers are not scavenged
564	flag_cvltr(it)=.FALSE.
565	ELSE IF(ANY(types_trac == 'inco')) THEN ! Add ThL
566	flag_cvltr(it)=.FALSE.
567	#ifdef CPP_StratAer
568	ELSE IF(ANY(types_trac == 'coag')) THEN
569	IF (convscav.and.aerosol(it)) THEN
570	flag_cvltr(it)=.TRUE.
571	ccntrAA(it) =ccntrAA_in
572	ccntrENV(it)=ccntrENV_in
573	coefcoli(it)=coefcoli_in
574	ELSE
575	flag_cvltr(it)=.FALSE.
576	ENDIF
577	#endif
578	ELSE IF(ANY(types_trac == 'lmdz')) THEN
579	IF (convscav.and.aerosol(it)) THEN
580	flag_cvltr(it)=.TRUE.
581	ccntrAA(it) =ccntrAA_in !--a modifier par JYG a lire depuis fichier
582	ccntrENV(it)=ccntrENV_in
583	coefcoli(it)=coefcoli_in
584	ELSE
585	flag_cvltr(it)=.FALSE.
586	ENDIF
587	END IF
588	ENDDO
589	!
590	ELSE ! iflag_con .ne. 3
591	flag_cvltr(:) = .FALSE.
592	ENDIF
593	!
594	! print out all tracer flags
595	!
596	WRITE(lunout,*) 'print out all tracer flags'
597	WRITE(lunout,*) 'type_trac =', type_trac
598	WRITE(lunout,*) 'config_inca =', config_inca
599	WRITE(lunout,*) 'iflag_con_trac =', iflag_con_trac
600	WRITE(lunout,*) 'iflag_con =', iflag_con
601	WRITE(lunout,*) 'convscav =', convscav
602	WRITE(lunout,*) 'iflag_lscav =', iflag_lscav
603	WRITE(lunout,*) 'aerosol =', aerosol
604	WRITE(lunout,*) 'iflag_the_trac =', iflag_the_trac
605	WRITE(lunout,*) 'iflag_thermals =', iflag_thermals
606	WRITE(lunout,*) 'iflag_vdf_trac =', iflag_vdf_trac
607	WRITE(lunout,*) 'pbl_flg =', pbl_flg
608	WRITE(lunout,*) 'lessivage =', lessivage
609	write(lunout,*) 'flag_cvltr = ', flag_cvltr
610
611	IF (lessivage .AND. (ANY(types_trac == 'inca') .OR. ANY(types_trac=='inco'))) THEN ! Mod ThL
612	CALL abort_physic('phytrac', 'lessivage=T config_inca=inca impossible',1)
613	! STOP
614	ENDIF
615	!
616	ENDIF ! of IF (debutphy)
617	!############################################ END INITIALIZATION #######
618
619	DO k=1,klev
620	DO i=1,klon
621	zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
622	END DO
623	END DO
624	!
625	IF (id_be .GT. 0) THEN
626	DO k=1,klev
627	DO i=1,klon
628	sourceBE(i,k)=srcbe(i,k) !RomP -> pour sortie histrac
629	END DO
630	END DO
631	ENDIF
632
633	!===============================================================================
634	! -- Do specific treatment according to chemestry model or local LMDZ tracers
635	!
636	!===============================================================================
637	IF(ANY(types_trac == 'inca')) THEN
638	! -- CHIMIE INCA config_inca = aero or chem --
639	! Appel fait en fin de phytrac pour avoir les emissions modifiees par
640	! la couche limite et la convection avant le calcul de la chimie
641
642	ELSE IF(ANY(types_trac == 'repr')) THEN
643	! -- CHIMIE REPROBUS --
644	CALL tracreprobus(pdtphys, gmtime, debutphy, julien, &
645	presnivs, xlat, xlon, pphis, pphi, &
646	t_seri, pplay, paprs, sh , &
647	tr_seri)
648
649	ELSE IF(ANY(types_trac == 'co2i')) THEN
650	! -- CO2 interactif --
651	! -- source is updated with FF and BB emissions
652	! -- and net fluxes from ocean and orchidee
653	! -- sign convention : positive into the atmosphere
654
655	CALL tracco2i(pdtphys, debutphy, &
656	xlat, xlon, pphis, pphi, &
657	t_seri, pplay, paprs, tr_seri, source)
658	ELSE IF(ANY(types_trac == 'inco')) THEN ! Add ThL
659	CALL tracco2i(pdtphys, debutphy, &
660	xlat, xlon, pphis, pphi, &
661	t_seri, pplay, paprs, tr_seri, source)
662
663	#ifdef CPP_StratAer
664	ELSE IF(ANY(types_trac == 'coag')) THEN
665	! --STRATOSPHERIC AER IN THE STRAT --
666	CALL traccoag(pdtphys, gmtime, debutphy, julien, &
667	presnivs, xlat, xlon, pphis, pphi, &
668	t_seri, pplay, paprs, sh, rh , &
669	tr_seri)
670	#endif
671	ELSE IF(ANY(types_trac == 'lmdz')) THEN
672	! -- Traitement des traceurs avec traclmdz
673	CALL traclmdz(nstep, julien, gmtime, pdtphys, t_seri, paprs, pplay, &
674	cdragh, coefh, yu1, yv1, ftsol, pctsrf, xlat, xlon,iflag_vdf_trac>=0,sh, &
675	rh, pphi, ustar, wstar, ale_bl, ale_wake, u10m, v10m, &
676	tr_seri, source, d_tr_cl,d_tr_dec, zmasse) !RomP
677	END IF
678	!======================================================================
679	! -- Calcul de l'effet de la convection --
680	!======================================================================
681
682	IF (iflag_con_trac==1) THEN
683
684	DO it=1, nbtr
685	IF ( conv_flg(it) == 0 ) CYCLE
686	IF (iflag_con.LT.2) THEN
687	!--pas de transport convectif
688	d_tr_cv(:,:,it)=0.
689
690	ELSE IF (iflag_con.EQ.2) THEN
691	!--ancien transport convectif de Tiedtke
692
693	CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
694	pplay, paprs, tr_seri(:,:,it), d_tr_cv(:,:,it))
695	ELSE
696	!--nouveau transport convectif de Emanuel
697
698	IF (flag_cvltr(it)) THEN
699	!--nouveau transport convectif de Emanuel avec lessivage convectif
700	!
701	!
702	ccntrAA_3d(:,:) =ccntrAA(it)
703	ccntrENV_3d(:,:)=ccntrENV(it)
704	coefcoli_3d(:,:)=coefcoli(it)
705
706	!--beware this interface is a bit weird because it is called for each tracer
707	!--with the full array tr_seri even if only item it is processed
708
709	print*,'CV SCAV ',it,ccntrAA(it),ccntrENV(it)
710
711	CALL cvltr_scav(pdtphys, da, phi,phi2,d1a,dam, mp,ep, &
712	sigd,sij,wght_cvfd,clw,elij,epmlmMm,eplaMm, &
713	pmflxr,pmflxs,evap,t_seri,wdtrainA,wdtrainM, &
714	paprs,it,tr_seri,upwd,dnwd,itop_con,ibas_con, &
715	ccntrAA_3d,ccntrENV_3d,coefcoli_3d, &
716	d_tr_cv,d_tr_trsp,d_tr_sscav,d_tr_sat,d_tr_uscav,qDi,qPr,&
717	qPa,qMel,qTrdi,dtrcvMA,Mint, &
718	zmfd1a,zmfphi2,zmfdam)
719
720
721	ELSE !---flag_cvltr(it).EQ.FALSE
722	!--nouveau transport convectif de Emanuel mais pas de lessivage convectif
723
724	!--beware this interface is a bit weird because it is called for each tracer
725	!--with the full array tr_seri even if only item it is processed
726	!
727	CALL cvltr_noscav(it,pdtphys, da, phi,mp,wght_cvfd,paprs,pplay, & !jyg
728	tr_seri,upwd,dnwd,d_tr_cv) !jyg
729
730	ENDIF
731
732	ENDIF !--iflag
733
734	!--on ajoute les tendances
735
736	DO k = 1, klev
737	DO i = 1, klon
738	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it)
739	END DO
740	END DO
741
742	CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'convection it = '//TRIM(int2str(it)))
743
744	END DO ! nbtr
745
746	#ifdef CPP_StratAer
747	IF (ANY(types_trac=='coag')) THEN
748	! initialize wet deposition flux of sulfur
749	budg_dep_wet_ocs(:)=0.0
750	budg_dep_wet_so2(:)=0.0
751	budg_dep_wet_h2so4(:)=0.0
752	budg_dep_wet_part(:)=0.0
753	! compute wet deposition flux of sulfur (sum over gases and particles)
754	! and convert to kg(S)/m2/s
755	DO i = 1, klon
756	DO k = 1, klev
757	DO it = 1, nbtr
758	!do not include SO2 because most of it comes trom the troposphere
759	IF (it==id_OCS_strat) THEN
760	budg_dep_wet_ocs(i)=budg_dep_wet_ocs(i)+d_tr_cv(i,k,it)*(mSatom/mOCSmol) &
761	& *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
762	ELSEIF (it==id_SO2_strat) THEN
763	budg_dep_wet_so2(i)=budg_dep_wet_so2(i)+d_tr_cv(i,k,it)*(mSatom/mSO2mol) &
764	& *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
765	ELSEIF (it==id_H2SO4_strat) THEN
766	budg_dep_wet_h2so4(i)=budg_dep_wet_h2so4(i)+d_tr_cv(i,k,it)*(mSatom/mH2SO4mol) &
767	& *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
768	ELSEIF (it.GT.nbtr_sulgas) THEN
769	budg_dep_wet_part(i)=budg_dep_wet_part(i)+d_tr_cv(i,k,it)*(mSatom/mH2SO4mol) &
770	& dens_aer_dry4./3.RPI(mdw(it-nbtr_sulgas)/2.)**3 &
771	& *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
772	ENDIF
773	ENDDO
774	ENDDO
775	ENDDO
776	ENDIF
777	#endif
778
779	ENDIF ! convection
780
781	!======================================================================
782	! -- Calcul de l'effet des thermiques --
783	!======================================================================
784
785	DO it=1,nbtr
786	DO k=1,klev
787	DO i=1,klon
788	d_tr_th(i,k,it)=0.
789	tr_seri(i,k,it)=MAX(tr_seri(i,k,it),0.)
790	! the next safeguard causes some problem for stratospheric aerosol tracers (particle number)
791	! and there is little justification for it so it is commented out (4 December 2017) by OB
792	! if reinstated please keep the ifndef CPP_StratAer
793	!#ifndef CPP_StratAer
794	! tr_seri(i,k,it)=MIN(tr_seri(i,k,it),1.e10)
795	!#endif
796	END DO
797	END DO
798	END DO
799
800	IF (iflag_thermals.GT.0.AND.iflag_the_trac>0) THEN
801
802	DO it=1, nbtr
803
804	CALL thermcell_dq(klon,klev,1,pdtphys,fm_therm,entr_therm, &
805	zmasse,tr_seri(1:klon,1:klev,it), &
806	d_tr_th(1:klon,1:klev,it),ztra_th,0 )
807
808	DO k=1,klev
809	DO i=1,klon
810	d_tr_th(i,k,it)=pdtphys*d_tr_th(i,k,it)
811	tr_seri(i,k,it)=MAX(tr_seri(i,k,it)+d_tr_th(i,k,it),0.)
812	END DO
813	END DO
814
815	END DO ! it
816
817	ENDIF ! Thermiques
818
819	!======================================================================
820	! -- Calcul de l'effet de la couche limite --
821	!======================================================================
822
823	IF (iflag_vdf_trac==1) THEN
824
825	! Injection during BL mixing
826	!
827	#ifdef CPP_StratAer
828	IF (ANY(types_trac=='coag')) THEN
829
830	! initialize dry deposition flux of sulfur
831	budg_dep_dry_ocs(:)=0.0
832	budg_dep_dry_so2(:)=0.0
833	budg_dep_dry_h2so4(:)=0.0
834	budg_dep_dry_part(:)=0.0
835
836	! compute dry deposition velocity as function of surface type (numbers
837	! from IPSL note 23, 2002)
838	v_dep_dry(:) = pctsrf(:,is_ter) * 2.5e-3 &
839	& + pctsrf(:,is_oce) * 0.5e-3 &
840	& + pctsrf(:,is_lic) * 2.5e-3 &
841	& + pctsrf(:,is_sic) * 2.5e-3
842
843	! compute surface dry deposition flux
844	zrho(:,1)=pplay(:,1)/t_seri(:,1)/RD
845
846	DO it=1, nbtr
847	source(:,it) = - v_dep_dry(:) * tr_seri(:,1,it) * zrho(:,1)
848	ENDDO
849
850	ENDIF
851	#endif
852
853	DO it=1, nbtr
854	!
855	IF( pbl_flg(it) /= 0 ) THEN
856	!
857	CALL cltrac(pdtphys, coefh,t_seri, &
858	tr_seri(:,:,it), source(:,it), &
859	paprs, pplay, delp, &
860	d_tr_cl(:,:,it),d_tr_dry(:,it),flux_tr_dry(:,it))
861	!
862	tr_seri(:,:,it)=tr_seri(:,:,it)+d_tr_cl(:,:,it)
863	!
864	#ifdef CPP_StratAer
865	IF (ANY(types_trac=='coag')) THEN
866	! compute dry deposition flux of sulfur (sum over gases and particles)
867	IF (it==id_OCS_strat) THEN
868	budg_dep_dry_ocs(:)=budg_dep_dry_ocs(:)-source(:,it)*(mSatom/mOCSmol)
869	ELSEIF (it==id_SO2_strat) THEN
870	budg_dep_dry_so2(:)=budg_dep_dry_so2(:)-source(:,it)*(mSatom/mSO2mol)
871	ELSEIF (it==id_H2SO4_strat) THEN
872	budg_dep_dry_h2so4(:)=budg_dep_dry_h2so4(:)-source(:,it)*(mSatom/mH2SO4mol)
873	ELSEIF (it.GT.nbtr_sulgas) THEN
874	budg_dep_dry_part(:)=budg_dep_dry_part(:)-source(:,it)(mSatom/mH2SO4mol)dens_aer_dry &
875	& 4./3.RPI(mdw(it-nbtr_sulgas)/2.)*3
876	ENDIF
877	ENDIF
878	#endif
879	!
880	ENDIF
881	!
882	ENDDO
883	!
884	ELSE IF (iflag_vdf_trac==0) THEN
885	!
886	! Injection of source in the first model layer
887	!
888	DO it=1,nbtr
889	d_tr_cl(:,1,it)=source(:,it)RG/delp(:,1)pdtphys
890	tr_seri(:,1,it)=tr_seri(:,1,it)+d_tr_cl(:,1,it)
891	ENDDO
892	d_tr_cl(:,2:klev,1:nbtr)=0.
893	!
894	ELSE IF (iflag_vdf_trac==-1) THEN
895	!
896	! Nothing happens
897	d_tr_cl=0.
898	!
899	ELSE
900	!
901	CALL abort_physic('iflag_vdf_trac', 'cas non prevu',1)
902	!
903	ENDIF ! couche limite
904
905	!======================================================================
906	! Calcul de l'effet de la precipitation grande echelle
907	! POUR INCA le lessivage est fait directement dans INCA
908	!======================================================================
909
910	IF (lessivage) THEN
911
912	ql_incloud_ref = 10.e-4
913	ql_incloud_ref = 5.e-4
914
915
916	! calcul du contenu en eau liquide au sein du nuage
917	ql_incl = ql_incloud_ref
918	! choix du lessivage
919	!
920	IF (iflag_lscav .EQ. 3 .OR. iflag_lscav .EQ. 4) THEN
921	! ******** Olivier Boucher version (3) possibly with modified ql_incl (4)
922	!
923	DO it = 1, nbtr
924
925	IF (aerosol(it)) THEN
926	! incloud scavenging and removal by large scale rain ! orig : ql_incl was replaced by 0.5e-3 kg/kg
927	! the value 0.5e-3 kg/kg is from Giorgi and Chameides (1986), JGR
928	! Liu (2001) proposed to use 1.5e-3 kg/kg
929
930	!jyg<
931	!! CALL lsc_scav(pdtphys,it,iflag_lscav,ql_incl,prfl,psfl,rneb,beta_fisrt, &
932	CALL lsc_scav(pdtphys,it,iflag_lscav,aerosol,ql_incl,prfl,psfl,rneb,beta_fisrt, &
933	!>jyg
934	beta_v1,pplay,paprs,t_seri,tr_seri,d_tr_insc, &
935	d_tr_bcscav,d_tr_evapls,qPrls)
936
937	!large scale scavenging tendency
938	DO k = 1, klev
939	DO i = 1, klon
940	d_tr_ls(i,k,it)=d_tr_insc(i,k,it)+d_tr_bcscav(i,k,it)+d_tr_evapls(i,k,it)
941	tr_seri(i,k,it)=tr_seri(i,k,it)+d_tr_ls(i,k,it)
942	ENDDO
943	ENDDO
944	CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'lsc scav it = '//TRIM(int2str(it)))
945	ENDIF
946
947	END DO !tr
948
949	#ifdef CPP_StratAer
950	IF (ANY(types_trac=='coag')) THEN
951	! compute wet deposition flux of sulfur (sum over gases and
952	! particles) and convert to kg(S)/m2/s
953	! adding contribution of d_tr_ls to d_tr_cv (above)
954	DO i = 1, klon
955	DO k = 1, klev
956	DO it = 1, nbtr
957	IF (it==id_OCS_strat) THEN
958	budg_dep_wet_ocs(i)=budg_dep_wet_ocs(i)+d_tr_ls(i,k,it)*(mSatom/mOCSmol) &
959	& *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
960	ELSEIF (it==id_SO2_strat) THEN
961	budg_dep_wet_so2(i)=budg_dep_wet_so2(i)+d_tr_ls(i,k,it)*(mSatom/mSO2mol) &
962	& *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
963	ELSEIF (it==id_H2SO4_strat) THEN
964	budg_dep_wet_h2so4(i)=budg_dep_wet_h2so4(i)+d_tr_ls(i,k,it)*(mSatom/mH2SO4mol) &
965	& *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
966	ELSEIF (it.GT.nbtr_sulgas) THEN
967	budg_dep_wet_part(i)=budg_dep_wet_part(i)+d_tr_ls(i,k,it)*(mSatom/mH2SO4mol) &
968	& dens_aer_dry4./3.RPI(mdw(it-nbtr_sulgas)/2.)**3 &
969	& *(paprs(i,k)-paprs(i,k+1))/RG/pdtphys
970	ENDIF
971	ENDDO
972	ENDDO
973	ENDDO
974	ENDIF
975	#endif
976
977	ELSE IF (iflag_lscav .EQ. 2) THEN ! frac_impa, frac_nucl
978	! ********* modified old version
979
980	d_tr_lessi_nucl(:,:,:) = 0.
981	d_tr_lessi_impa(:,:,:) = 0.
982	flestottr(:,:,:) = 0.
983	! Tendance des aerosols nuclees et impactes
984	DO it = 1, nbtr
985	IF (aerosol(it)) THEN
986	his_dh(:)=0.
987	DO k = 1, klev
988	DO i = 1, klon
989	!PhH
990	zrho(i,k)=pplay(i,k)/t_seri(i,k)/RD
991	zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho(i,k)/RG
992	!
993	ENDDO
994	ENDDO
995
996	DO k=klev-1, 1, -1
997	DO i=1, klon
998	! d_tr_ls(i,k,it)=tr_seri(i,k,it)(frac_impa(i,k)frac_nucl(i,k)-1.)
999	dx=d_tr_ls(i,k,it)
1000	his_dh(i)=his_dh(i)-dxzrho(i,k)zdz(i,k)/pdtphys ! kg/m2/s
1001	evaplsc = prfl(i,k) - prfl(i,k+1) + psfl(i,k) - psfl(i,k+1)
1002	! Evaporation Partielle -> Liberation Partielle 0.5*evap
1003	IF ( evaplsc .LT.0..and.abs(prfl(i,k+1)+psfl(i,k+1)).gt.1.e-10) THEN
1004	evaplsc = (-evaplsc)/(prfl(i,k+1)+psfl(i,k+1))
1005	! evaplsc est donc positif, his_dh(i) est positif
1006	!--------------
1007	d_tr_evapls(i,k,it)=0.5evaplsc(d_tr_lessi_nucl(i,k+1,it) &
1008	+d_tr_lessi_impa(i,k+1,it))
1009	!------------- d_tr_evapls(i,k,it)=-0.5evaplsc(d_tr_lsc(i,k+1,it))
1010	beta=0.5*evaplsc
1011	if ((prfl(i,k)+psfl(i,k)).lt.1.e-10) THEN
1012	beta=1.0*evaplsc
1013	endif
1014	dx=betahis_dh(i)/zrho(i,k)/zdz(i,k)pdtphys
1015	his_dh(i)=(1.-beta)*his_dh(i) ! tracer from
1016	d_tr_evapls(i,k,it)=dx
1017	ENDIF
1018	d_tr_ls(i,k,it)=tr_seri(i,k,it)(frac_impa(i,k)frac_nucl(i,k)-1.) &
1019	+d_tr_evapls(i,k,it)
1020
1021	!--------------
1022	d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) + &
1023	( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
1024	d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) + &
1025	( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
1026	!
1027	! Flux lessivage total
1028	flestottr(i,k,it) = flestottr(i,k,it) - &
1029	( d_tr_lessi_nucl(i,k,it) + &
1030	d_tr_lessi_impa(i,k,it) ) * &
1031	( paprs(i,k)-paprs(i,k+1) ) / &
1032	(RG * pdtphys)
1033	!! Mise a jour des traceurs due a l'impaction,nucleation
1034	! tr_seri(i,k,it)=tr_seri(i,k,it)frac_impa(i,k)frac_nucl(i,k)
1035	!! calcul de la tendance liee au lessivage stratiforme
1036	! d_tr_ls(i,k,it)=tr_seri(i,k,it)*&
1037	! (1.-1./(frac_impa(i,k)*frac_nucl(i,k)))
1038	!--------------
1039	ENDDO
1040	ENDDO
1041	ENDIF
1042	ENDDO
1043	! ********* end modified old version
1044
1045	ELSE IF (iflag_lscav .EQ. 1) THEN ! frac_impa, frac_nucl
1046	! ********* old version
1047
1048	d_tr_lessi_nucl(:,:,:) = 0.
1049	d_tr_lessi_impa(:,:,:) = 0.
1050	flestottr(:,:,:) = 0.
1051	!=========================
1052	! LESSIVAGE LARGE SCALE :
1053	!=========================
1054
1055	! Tendance des aerosols nuclees et impactes
1056	! -----------------------------------------
1057	DO it = 1, nbtr
1058	IF (aerosol(it)) THEN
1059	DO k = 1, klev
1060	DO i = 1, klon
1061	d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) + &
1062	( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
1063	d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) + &
1064	( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
1065
1066	!
1067	! Flux lessivage total
1068	! ------------------------------------------------------------
1069	flestottr(i,k,it) = flestottr(i,k,it) - &
1070	( d_tr_lessi_nucl(i,k,it) + &
1071	d_tr_lessi_impa(i,k,it) ) * &
1072	( paprs(i,k)-paprs(i,k+1) ) / &
1073	(RG * pdtphys)
1074	!
1075	! Mise a jour des traceurs due a l'impaction,nucleation
1076	! ----------------------------------------------------------------------
1077	tr_seri(i,k,it)=tr_seri(i,k,it)frac_impa(i,k)frac_nucl(i,k)
1078	ENDDO
1079	ENDDO
1080	ENDIF
1081	ENDDO
1082
1083	! ********* end old version
1084	ENDIF ! iflag_lscav . EQ. 1, 2, 3 or 4
1085	!
1086	ENDIF ! lessivage
1087
1088
1089	! -- CHIMIE INCA config_inca = aero or chem --
1090	IF (ANY(types_trac == 'inca') .OR. ANY(types_trac == 'inco')) THEN ! ModThL
1091
1092	CALL tracinca(&
1093	nstep, julien, gmtime, lafin, &
1094	pdtphys, t_seri, paprs, pplay, &
1095	pmfu, upwd, ftsol, pctsrf, pphis, &
1096	pphi, albsol, sh, ch, rh, &
1097	cldfra, rneb, diafra, cldliq, &
1098	itop_con, ibas_con, pmflxr, pmflxs, &
1099	prfl, psfl, aerosol_couple, flxmass_w, &
1100	tau_aero, piz_aero, cg_aero, ccm, &
1101	rfname, &
1102	tr_seri(:,:,1+nqCO2:nbtr), source(:,1+nqCO2:nbtr)) ! ModThL
1103	ENDIF
1104
1105	END SUBROUTINE phytrac
1106
1107	END MODULE

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format