Context Navigation

phytrac_mod.F90 @ 2171

Last change on this file since 2171 was 2171, checked in by acozic, 10 years ago
There are some commits that we must not do just before holiday .... so be back to rev 2168
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
File size: 36.0 KB

Line
1	!$Id$
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,d_tr_cl,d_tr_dry,flux_tr_dry,d_tr_dec,d_tr_cv)
52
53
54	CONTAINS
55
56	SUBROUTINE phytrac( &
57	nstep, julien, gmtime, debutphy, &
58	lafin, pdtphys, u, v, t_seri, &
59	paprs, pplay, pmfu, pmfd, &
60	pen_u, pde_u, pen_d, pde_d, &
61	cdragh, coefh, fm_therm, entr_therm, &
62	yu1, yv1, ftsol, pctsrf, &
63	ustar, u10m, v10m, &
64	wstar, ale_bl, ale_wake, &
65	xlat, xlon, &
66	frac_impa,frac_nucl,beta_fisrt,beta_v1, &
67	presnivs, pphis, pphi, albsol, &
68	sh, rh, cldfra, rneb, &
69	diafra, cldliq, itop_con, ibas_con, &
70	pmflxr, pmflxs, prfl, psfl, &
71	da, phi, mp, upwd, &
72	phi2, d1a, dam, sij, wght_cvfd, & ! RomP +RL
73	wdtrainA, wdtrainM, sigd, clw, elij, & ! RomP
74	evap, ep, epmlmMm, eplaMm, & ! RomP
75	dnwd, aerosol_couple, flxmass_w, &
76	tau_aero, piz_aero, cg_aero, ccm, &
77	rfname, &
78	d_tr_dyn, & ! RomP
79	tr_seri)
80	!
81	!======================================================================
82	! Auteur(s) FH
83	! Objet: Moniteur general des tendances traceurs
84	! Modification R. Pilon 01 janvier 2012 transport+scavenging KE scheme : cvltr
85	! Modification R. Pilon 10 octobre 2012 large scale scavenging incloud_scav + bc_scav
86	!======================================================================
87
88	USE ioipsl
89	USE phys_cal_mod, only : hour
90	USE dimphy
91	USE infotrac
92	USE mod_grid_phy_lmdz
93	USE mod_phys_lmdz_para
94	USE comgeomphy
95	USE iophy
96	USE traclmdz_mod
97	USE tracinca_mod
98	USE tracreprobus_mod
99	USE control_mod
100	USE indice_sol_mod
101
102	USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root
103
104	IMPLICIT NONE
105
106	INCLUDE "YOMCST.h"
107	INCLUDE "dimensions.h"
108	INCLUDE "clesphys.h"
109	INCLUDE "temps.h"
110	INCLUDE "paramet.h"
111	INCLUDE "thermcell.h"
112	INCLUDE "iniprint.h"
113	!==========================================================================
114	! -- ARGUMENT DESCRIPTION --
115	!==========================================================================
116
117	! Input arguments
118	!----------------
119	!Configuration grille,temps:
120	INTEGER,INTENT(IN) :: nstep ! Appel physique
121	INTEGER,INTENT(IN) :: julien ! Jour julien
122	REAL,INTENT(IN) :: gmtime ! Heure courante
123	REAL,INTENT(IN) :: pdtphys ! Pas d'integration pour la physique (seconde)
124	LOGICAL,INTENT(IN) :: debutphy ! le flag de l'initialisation de la physique
125	LOGICAL,INTENT(IN) :: lafin ! le flag de la fin de la physique
126
127	REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point
128	REAL,DIMENSION(klon),INTENT(IN) :: xlon ! longitudes pour chaque point
129	!
130	!Physique:
131	!--------
132	REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature
133	REAL,DIMENSION(klon,klev),INTENT(IN) :: u ! variable not used
134	REAL,DIMENSION(klon,klev),INTENT(IN) :: v ! variable not used
135	REAL,DIMENSION(klon,klev),INTENT(IN) :: sh ! humidite specifique
136	REAL,DIMENSION(klon,klev),INTENT(IN) :: rh ! humidite relative
137	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa)
138	REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa)
139	REAL,DIMENSION(klon,klev),INTENT(IN) :: pphi ! geopotentiel
140	REAL,DIMENSION(klon),INTENT(IN) :: pphis
141	REAL,DIMENSION(klev),INTENT(IN) :: presnivs
142	REAL,DIMENSION(klon,klev),INTENT(IN) :: cldliq ! eau liquide nuageuse
143	REAL,DIMENSION(klon,klev),INTENT(IN) :: cldfra ! fraction nuageuse (tous les nuages)
144	REAL,DIMENSION(klon,klev),INTENT(IN) :: diafra ! fraction nuageuse (convection ou stratus artificiels)
145	REAL,DIMENSION(klon,klev),INTENT(IN) :: rneb ! fraction nuageuse (grande echelle)
146	!
147	REAL :: ql_incl ! contenu en eau liquide nuageuse dans le nuage ! ql_incl=oliq/rneb
148	REAL,DIMENSION(klon,klev),INTENT(IN) :: beta_fisrt ! taux de conversion de l'eau cond (de fisrtilp)
149	REAL,DIMENSION(klon,klev),INTENT(out) :: beta_v1 ! -- (originale version)
150
151	!
152	INTEGER,DIMENSION(klon),INTENT(IN) :: itop_con
153	INTEGER,DIMENSION(klon),INTENT(IN) :: ibas_con
154	REAL,DIMENSION(klon),INTENT(IN) :: albsol ! albedo surface
155	!
156	!Dynamique
157	!--------
158	REAL,DIMENSION(klon,klev,nbtr),INTENT(IN) :: d_tr_dyn
159	!
160	!Convection:
161	!----------
162	REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfu ! flux de masse dans le panache montant
163	REAL,DIMENSION(klon,klev),INTENT(IN) :: pmfd ! flux de masse dans le panache descendant
164	REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_u ! flux entraine dans le panache montant
165	REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_u ! flux detraine dans le panache montant
166	REAL,DIMENSION(klon,klev),INTENT(IN) :: pen_d ! flux entraine dans le panache descendant
167	REAL,DIMENSION(klon,klev),INTENT(IN) :: pde_d ! flux detraine dans le panache descendant
168
169	!...Tiedke
170	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: pmflxr, pmflxs ! Flux precipitant de pluie, neige aux interfaces [convection]
171	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: prfl, psfl ! Flux precipitant de pluie, neige aux interfaces [large-scale]
172
173	LOGICAL,INTENT(IN) :: aerosol_couple
174	REAL,DIMENSION(klon,klev),INTENT(IN) :: flxmass_w
175	REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: tau_aero
176	REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: piz_aero
177	REAL,DIMENSION(klon,klev,9,2),INTENT(IN) :: cg_aero
178	CHARACTER(len=4),DIMENSION(9),INTENT(IN) :: rfname
179	REAL,DIMENSION(klon,klev,2),INTENT(IN) :: ccm
180	!... K.Emanuel
181	REAL,DIMENSION(klon,klev),INTENT(IN) :: da
182	REAL,DIMENSION(klon,klev,klev),INTENT(IN):: phi
183	! RomP >>>
184	REAL,DIMENSION(klon,klev),INTENT(IN) :: d1a,dam
185	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: phi2
186	!
187	REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtrainA
188	REAL,DIMENSION(klon,klev),INTENT(IN) :: wdtrainM
189	REAL,DIMENSION(klon),INTENT(IN) :: sigd
190	! ---- RomP flux entraine, detraine et precipitant kerry Emanuel
191	REAL,DIMENSION(klon,klev),INTENT(IN) :: evap
192	REAL,DIMENSION(klon,klev),INTENT(IN) :: ep
193	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: sij
194	REAL,DIMENSION(klon,klev),INTENT(IN) :: wght_cvfd !RL
195	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: elij
196	REAL,DIMENSION(klon,klev,klev),INTENT(IN) :: epmlmMm
197	REAL,DIMENSION(klon,klev),INTENT(IN) :: eplaMm
198	REAL,DIMENSION(klon,klev),INTENT(IN) :: clw
199	! RomP <<<
200
201	!
202	REAL,DIMENSION(klon,klev),INTENT(IN) :: mp
203	REAL,DIMENSION(klon,klev),INTENT(IN) :: upwd ! saturated updraft mass flux
204	REAL,DIMENSION(klon,klev),INTENT(IN) :: dnwd ! saturated downdraft mass flux
205	!
206	!Thermiques:
207	!----------
208	REAL,DIMENSION(klon,klev+1),INTENT(IN) :: fm_therm
209	REAL,DIMENSION(klon,klev),INTENT(IN) :: entr_therm
210	!
211	!Couche limite:
212	!--------------
213	!
214	!
215	REAL,DIMENSION(:),INTENT(IN) :: cdragh ! (klon) coeff drag pour T et Q
216	REAL,DIMENSION(:,:),INTENT(IN) :: coefh ! (klon,klev) coeff melange CL (m**2/s)
217	REAL,DIMENSION(:),INTENT(IN) :: ustar,u10m,v10m ! (klon) u* & vent a 10m (m/s)
218	REAL,DIMENSION(:),INTENT(IN) :: wstar,ale_bl,ale_wake ! (klon) w* and Avail. Lifting Ener.
219	REAL,DIMENSION(:),INTENT(IN) :: yu1 ! (klon) vents au premier niveau
220	REAL,DIMENSION(:),INTENT(IN) :: yv1 ! (klon) vents au premier niveau
221
222	!
223	!Lessivage:
224	!----------
225	!
226	REAL, DIMENSION(:), ALLOCATABLE, SAVE :: ccntrAA
227	REAL, DIMENSION(:), ALLOCATABLE, SAVE :: ccntrENV
228	REAL, DIMENSION(:), ALLOCATABLE, SAVE :: coefcoli
229	LOGICAL, DIMENSION(:), ALLOCATABLE, SAVE :: flag_cvltr
230	!$OMP THREADPRIVATE(ccntrAA,ccntrENV,coefcoli,flag_cvltr)
231	REAL, DIMENSION(klon,klev) :: ccntrAA_3d
232	REAL, DIMENSION(klon,klev) :: ccntrENV_3d
233	REAL, DIMENSION(klon,klev) :: coefcoli_3d
234	!
235	! pour le ON-LINE
236	!
237	REAL,DIMENSION(klon,klev),INTENT(IN) :: frac_impa ! fraction d'aerosols non impactes
238	REAL,DIMENSION(klon,klev),INTENT(IN) :: frac_nucl ! fraction d'aerosols non nuclees
239
240	! Arguments necessaires pour les sources et puits de traceur:
241	REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: ftsol ! Temperature du sol (surf)(Kelvin)
242	REAL,DIMENSION(klon,nbsrf),INTENT(IN) :: pctsrf ! Pourcentage de sol (nature du sol)
243
244
245	! Output argument
246	!----------------
247	REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA]
248	REAL,DIMENSION(klon,klev) :: sourceBE
249	!=======================================================================================
250	! -- LOCAL VARIABLES --
251	!=======================================================================================
252
253	INTEGER :: i, k, it
254	INTEGER :: nsplit
255
256	!Sources et Reservoirs de traceurs (ex:Radon):
257	!--------------------------------------------
258	!
259	REAL,DIMENSION(:,:),ALLOCATABLE,SAVE :: source ! a voir lorsque le flux de surface est prescrit
260	!$OMP THREADPRIVATE(source)
261
262	!
263	!Entrees/Sorties: (cf ini_histrac.h et write_histrac.h)
264	!---------------
265	INTEGER :: iiq, ierr
266	INTEGER :: nhori, nvert
267	REAL :: zsto, zout, zjulian
268	INTEGER,SAVE :: nid_tra ! pointe vers le fichier histrac.nc
269	!$OMP THREADPRIVATE(nid_tra)
270	REAL,DIMENSION(klon) :: zx_tmp_fi2d ! variable temporaire grille physique
271	INTEGER :: itau_w ! pas de temps ecriture = nstep + itau_phy
272	LOGICAL,PARAMETER :: ok_sync=.TRUE.
273
274	!
275	! Nature du traceur
276	!------------------
277	LOGICAL,DIMENSION(:),ALLOCATABLE,SAVE :: aerosol ! aerosol(it) = true => aerosol => lessivage
278	!$OMP THREADPRIVATE(aerosol) ! aerosol(it) = false => gaz
279	REAL,DIMENSION(klon,klev) :: delp ! epaisseur de couche (Pa)
280	!
281	! Tendances de traceurs (Td) et flux de traceurs:
282	!------------------------
283	REAL,DIMENSION(klon,klev) :: d_tr ! Td dans l'atmosphere
284	REAL,DIMENSION(klon,klev) :: Mint
285	REAL,DIMENSION(klon,klev,nbtr) :: zmfd1a
286	REAL,DIMENSION(klon,klev,nbtr) :: zmfdam
287	REAL,DIMENSION(klon,klev,nbtr) :: zmfphi2
288	! Physique
289	!----------
290	REAL,DIMENSION(klon,klev,nbtr) :: flestottr ! flux de lessivage dans chaque couche
291	REAL,DIMENSION(klon,klev) :: zmasse ! densité atmosphérique Kg/m2
292	REAL,DIMENSION(klon,klev) :: ztra_th
293	!PhH
294	REAL,DIMENSION(klon,klev) :: zrho
295	REAL,DIMENSION(klon,klev) :: zdz
296	REAL :: evaplsc,dx,beta ! variable pour lessivage Genthon
297	REAL,DIMENSION(klon) :: his_dh ! ---
298	! in-cloud scav variables
299	REAL :: ql_incloud_ref ! ref value of in-cloud condensed water content
300
301	!Controles:
302	!---------
303	INTEGER,SAVE :: iflag_vdf_trac,iflag_con_trac,iflag_the_trac
304	INTEGER,SAVE :: iflag_con_trac_omp, iflag_vdf_trac_omp,iflag_the_trac_omp
305	!$OMP THREADPRIVATE(iflag_vdf_trac,iflag_con_trac,iflag_the_trac)
306
307	LOGICAL,SAVE :: lessivage
308	!$OMP THREADPRIVATE(lessivage)
309
310	CHARACTER(len=8),DIMENSION(nbtr) :: solsym
311	!RomP >>>
312	INTEGER,SAVE :: iflag_lscav_omp,iflag_lscav
313	LOGICAL,SAVE :: convscav_omp,convscav
314	!$OMP THREADPRIVATE(iflag_lscav)
315	!$OMP THREADPRIVATE(convscav)
316	!RomP <<<
317	!######################################################################
318	! -- INITIALIZATION --
319	!######################################################################
320	IF (debutphy) THEN
321	ALLOCATE(d_tr_cl(klon,klev,nbtr),d_tr_dry(klon,nbtr))
322	ALLOCATE(flux_tr_dry(klon,nbtr),d_tr_dec(klon,klev,nbtr),d_tr_cv(klon,klev,nbtr))
323	ALLOCATE(d_tr_insc(klon,klev,nbtr),d_tr_bcscav(klon,klev,nbtr))
324	ALLOCATE(d_tr_evapls(klon,klev,nbtr),d_tr_ls(klon,klev,nbtr))
325	ALLOCATE(qPrls(klon,nbtr),d_tr_trsp(klon,klev,nbtr))
326	ALLOCATE(d_tr_sscav(klon,klev,nbtr),d_tr_sat(klon,klev,nbtr))
327	ALLOCATE(d_tr_uscav(klon,klev,nbtr),qPr(klon,klev,nbtr),qDi(klon,klev,nbtr))
328	ALLOCATE(qPa(klon,klev,nbtr),qMel(klon,klev,nbtr))
329	ALLOCATE(qTrdi(klon,klev,nbtr),dtrcvMA(klon,klev,nbtr))
330	ALLOCATE(d_tr_th(klon,klev,nbtr))
331	ALLOCATE(d_tr_lessi_impa(klon,klev,nbtr),d_tr_lessi_nucl(klon,klev,nbtr))
332	ENDIF
333
334	DO k=1,klev
335	DO i=1,klon
336	sourceBE(i,k)=0.
337	Mint(i,k)=0.
338	zrho(i,k)=0.
339	zdz(i,k)=0.
340	END DO
341	END DO
342
343	DO it=1, nbtr
344	DO k=1,klev
345	DO i=1,klon
346	d_tr_insc(i,k,it)=0.
347	d_tr_bcscav(i,k,it)=0.
348	d_tr_evapls(i,k,it)=0.
349	d_tr_ls(i,k,it)=0.
350	d_tr_cv(i,k,it)=0.
351	d_tr_cl(i,k,it)=0.
352	d_tr_trsp(i,k,it)=0.
353	d_tr_sscav(i,k,it)=0.
354	d_tr_sat(i,k,it)=0.
355	d_tr_uscav(i,k,it)=0.
356	d_tr_lessi_impa(i,k,it)=0.
357	d_tr_lessi_nucl(i,k,it)=0.
358	qDi(i,k,it)=0.
359	qPr(i,k,it)=0.
360	qPa(i,k,it)=0.
361	qMel(i,k,it)=0.
362	qTrdi(i,k,it)=0.
363	dtrcvMA(i,k,it)=0.
364	zmfd1a(i,k,it)=0.
365	zmfdam(i,k,it)=0.
366	zmfphi2(i,k,it)=0.
367	END DO
368	END DO
369	END DO
370
371	DO k = 1, klev
372	DO i = 1, klon
373	delp(i,k) = paprs(i,k)-paprs(i,k+1)
374	END DO
375	END DO
376
377	IF (debutphy) THEN
378	!!jyg
379	!$OMP BARRIER
380	ecrit_tra=86400. ! frequence de stokage en dur
381	! obsolete car remplace par des ecritures dans phys_output_write
382	!RomP >>>
383	!
384	!Config Key = convscav
385	!Config Desc = Convective scavenging switch: 0=off, 1=on.
386	!Config Def = .false.
387	!Config Help =
388	!
389	!$OMP MASTER
390	convscav_omp=.false.
391	call getin('convscav', convscav_omp)
392	iflag_vdf_trac_omp=1
393	call getin('iflag_vdf_trac', iflag_vdf_trac_omp)
394	iflag_con_trac_omp=1
395	call getin('iflag_con_trac', iflag_con_trac_omp)
396	iflag_the_trac_omp=1
397	call getin('iflag_the_trac', iflag_the_trac_omp)
398	!$OMP END MASTER
399	!$OMP BARRIER
400	convscav=convscav_omp
401	iflag_vdf_trac=iflag_vdf_trac_omp
402	iflag_con_trac=iflag_con_trac_omp
403	iflag_the_trac=iflag_the_trac_omp
404	write(lunout,*) 'phytrac passage dans routine conv avec lessivage', convscav
405	!
406	!Config Key = iflag_lscav
407	!Config Desc = Large scale scavenging parametrization: 0=none, 1=old(Genthon92),
408	! 2=1+PHeinrich, 3=Reddy_Boucher2004, 4=3+RPilon.
409	!Config Def = 1
410	!Config Help =
411	!
412	!$OMP MASTER
413	iflag_lscav_omp=1
414	call getin('iflag_lscav', iflag_lscav_omp)
415	!$OMP END MASTER
416	!$OMP BARRIER
417	iflag_lscav=iflag_lscav_omp
418	!
419	SELECT CASE(iflag_lscav)
420	CASE(0)
421	WRITE(lunout,*) 'Large scale scavenging: none'
422	CASE(1)
423	WRITE(lunout,*) 'Large scale scavenging: C. Genthon, Tellus(1992), 44B, 371-389'
424	CASE(2)
425	WRITE(lunout,*) 'Large scale scavenging: C. Genthon, modified P. Heinrich'
426	CASE(3)
427	WRITE(lunout,*) 'Large scale scavenging: M. Shekkar Reddy and O. Boucher, JGR(2004), 109, D14202'
428	CASE(4)
429	WRITE(lunout,*) 'Large scale scavenging: Reddy and Boucher, modified R. Pilon'
430	END SELECT
431	!RomP <<<
432	WRITE(,) 'FIRST TIME IN PHYTRAC : pdtphys(sec) = ',pdtphys,'ecrit_tra (sec) = ',ecrit_tra
433	ALLOCATE( source(klon,nbtr), stat=ierr)
434	IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 1',1)
435
436	ALLOCATE( aerosol(nbtr), stat=ierr)
437	IF (ierr /= 0) CALL abort_gcm('phytrac', 'pb in allocation 2',1)
438
439
440	! Initialize module for specific tracers
441	SELECT CASE(type_trac)
442	CASE('lmdz')
443	CALL traclmdz_init(pctsrf,xlat,xlon,ftsol,tr_seri,t_seri,pplay,sh,pdtphys,aerosol,lessivage)
444	CASE('inca')
445	source(:,:)=0.
446	CALL tracinca_init(aerosol,lessivage)
447	CASE('repr')
448	source(:,:)=0.
449	END SELECT
450
451	!
452	!--initialising coefficients for scavenging in the case of NP
453	!
454	ALLOCATE(flag_cvltr(nbtr))
455	IF (iflag_con.EQ.3) THEN
456	!
457	ALLOCATE(ccntrAA(nbtr))
458	ALLOCATE(ccntrENV(nbtr))
459	ALLOCATE(coefcoli(nbtr))
460	!
461	DO it=1, nbtr
462	SELECT CASE(type_trac)
463	CASE('lmdz')
464	IF (convscav.and.aerosol(it)) THEN
465	flag_cvltr(it)=.true.
466	ccntrAA(it) =1.0 !--a modifier par JYG a lire depuis fichier
467	ccntrENV(it)=1.0
468	coefcoli(it)=0.001
469	ELSE
470	flag_cvltr(it)=.false.
471	ENDIF
472
473	CASE('inca')
474	! IF ((it.EQ.id_Rn222) .OR. ((it.GE.id_SO2) .AND. (it.LE.id_NH3)) ) THEN
475	! !--gas-phase species
476	! flag_cvltr(it)=.false.
477	!
478	! ELSEIF ( (it.GE.id_CIDUSTM) .AND. (it.LE.id_AIN) ) THEN
479	! !--insoluble aerosol species
480	! flag_cvltr(it)=.true.
481	! ccntrAA(it)=0.7
482	! ccntrENV(it)=0.7
483	! coefcoli(it)=0.001
484	! ELSEIF ( (it.EQ.id_Pb210) .OR. ((it.GE.id_CSSSM) .AND. (it.LE.id_SSN))) THEN
485	! !--soluble aerosol species
486	! flag_cvltr(it)=.true.
487	! ccntrAA(it)=0.9
488	! ccntrENV(it)=0.9
489	! coefcoli(it)=0.001
490	! ELSE
491	! WRITE(lunout,*) 'pb it=', it
492	! CALL abort_gcm('phytrac','pb it scavenging',1)
493	! ENDIF
494	!--test OB
495	!--for now we do not scavenge in cvltr
496	flag_cvltr(it)=.false.
497	END SELECT
498	ENDDO
499	!
500	ELSE ! iflag_con .ne. 3
501	flag_cvltr(:) = .false.
502	ENDIF
503	!
504	! Initialize diagnostic output
505	! ----------------------------
506	#ifdef CPP_IOIPSL
507	! INCLUDE "ini_histrac.h"
508	#endif
509	!
510	! print out all tracer flags
511	!
512	WRITE(lunout,*) 'print out all tracer flags'
513	WRITE(lunout,*) 'type_trac =', type_trac
514	WRITE(lunout,*) 'config_inca =', config_inca
515	WRITE(lunout,*) 'iflag_con_trac =', iflag_con_trac
516	WRITE(lunout,*) 'iflag_con =', iflag_con
517	WRITE(lunout,*) 'convscav =', convscav
518	WRITE(lunout,*) 'iflag_lscav =', iflag_lscav
519	WRITE(lunout,*) 'aerosol =', aerosol
520	WRITE(lunout,*) 'iflag_the_trac =', iflag_the_trac
521	WRITE(lunout,*) 'iflag_thermals =', iflag_thermals
522	WRITE(lunout,*) 'iflag_vdf_trac =', iflag_vdf_trac
523	WRITE(lunout,*) 'pbl_flg =', pbl_flg
524	WRITE(lunout,*) 'lessivage =', lessivage
525	write(lunout,*) 'flag_cvltr = ', flag_cvltr
526
527	IF (lessivage.AND.config_inca.EQ.'inca') THEN
528	CALL abort_gcm('phytrac', 'lessivage=T config_inca=inca impossible',1)
529	STOP
530	ENDIF
531	!
532	END IF ! of IF (debutphy)
533	!############################################ END INITIALIZATION #######
534
535	DO k=1,klev
536	DO i=1,klon
537	zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
538	END DO
539	END DO
540	!
541	IF (id_be .GT. 0) THEN
542	DO k=1,klev
543	DO i=1,klon
544	sourceBE(i,k)=srcbe(i,k) !RomP -> pour sortie histrac
545	END DO
546	END DO
547	ENDIF
548
549	!===============================================================================
550	! -- Do specific treatment according to chemestry model or local LMDZ tracers
551	!
552	!===============================================================================
553	SELECT CASE(type_trac)
554	CASE('lmdz')
555	! -- Traitement des traceurs avec traclmdz
556	CALL traclmdz(nstep, julien, gmtime, pdtphys, t_seri, paprs, pplay, &
557	cdragh, coefh, yu1, yv1, ftsol, pctsrf, xlat, xlon,iflag_vdf_trac>=0,sh, &
558	rh, pphi, ustar, wstar, ale_bl, ale_wake, u10m, v10m, &
559	tr_seri, source, solsym, d_tr_cl,d_tr_dec, zmasse) !RomP
560
561	CASE('inca')
562	! -- CHIMIE INCA config_inca = aero or chem --
563
564	CALL tracinca(&
565	nstep, julien, gmtime, lafin, &
566	pdtphys, t_seri, paprs, pplay, &
567	pmfu, upwd, ftsol, pctsrf, pphis, &
568	pphi, albsol, sh, rh, &
569	cldfra, rneb, diafra, cldliq, &
570	itop_con, ibas_con, pmflxr, pmflxs, &
571	prfl, psfl, aerosol_couple, flxmass_w, &
572	tau_aero, piz_aero, cg_aero, ccm, &
573	rfname, &
574	tr_seri, source, solsym)
575
576	CASE('repr')
577	! -- CHIMIE REPROBUS --
578
579	CALL tracreprobus(pdtphys, gmtime, debutphy, julien, &
580	presnivs, xlat, xlon, pphis, pphi, &
581	t_seri, pplay, paprs, sh , &
582	tr_seri, solsym)
583
584	END SELECT
585	!======================================================================
586	! -- Calcul de l'effet de la convection --
587	!======================================================================
588
589	IF (iflag_con_trac==1) THEN
590
591	DO it=1, nbtr
592	IF ( conv_flg(it) == 0 ) CYCLE
593	IF (iflag_con.LT.2) THEN
594	!--pas de transport convectif
595
596	d_tr_cv(:,:,it)=0.
597	ELSE IF (iflag_con.EQ.2) THEN
598	!--ancien transport convectif de Tiedtke
599
600	CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d, &
601	pplay, paprs, tr_seri(:,:,it), d_tr_cv(:,:,it))
602	ELSE
603	!--nouveau transport convectif de Emanuel
604
605	IF (flag_cvltr(it)) THEN
606	!--nouveau transport convectif de Emanuel avec lessivage convectif
607	!
608	!
609	ccntrAA_3d(:,:) =ccntrAA(it)
610	ccntrENV_3d(:,:)=ccntrENV(it)
611	coefcoli_3d(:,:)=coefcoli(it)
612
613	!--beware this interface is a bit weird because it is called for each tracer
614	!--with the full array tr_seri even if only item it is processed
615
616	CALL cvltr_scav(pdtphys, da, phi,phi2,d1a,dam, mp,ep, &
617	sigd,sij,wght_cvfd,clw,elij,epmlmMm,eplaMm, &
618	pmflxr,pmflxs,evap,t_seri,wdtrainA,wdtrainM, &
619	paprs,it,tr_seri,upwd,dnwd,itop_con,ibas_con, &
620	ccntrAA_3d,ccntrENV_3d,coefcoli_3d, &
621	d_tr_cv,d_tr_trsp,d_tr_sscav,d_tr_sat,d_tr_uscav,qDi,qPr,&
622	qPa,qMel,qTrdi,dtrcvMA,Mint, &
623	zmfd1a,zmfphi2,zmfdam)
624
625
626	ELSE !---flag_cvltr(it).EQ.FALSE
627	!--nouveau transport convectif de Emanuel mais pas de lessivage convectif
628
629	!--beware this interface is a bit weird because it is called for each tracer
630	!--with the full array tr_seri even if only item it is processed
631	!
632	CALL cvltr_noscav(it,pdtphys, da, phi,mp,wght_cvfd,paprs,pplay, & !jyg
633	tr_seri,upwd,dnwd,d_tr_cv) !jyg
634
635	ENDIF
636
637	ENDIF !--iflag
638
639	!--on ajoute les tendances
640
641	DO k = 1, klev
642	DO i = 1, klon
643	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it)
644	END DO
645	END DO
646
647	CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'convection it = '//solsym(it))
648
649	END DO ! nbtr
650
651	END IF ! convection
652
653	!======================================================================
654	! -- Calcul de l'effet des thermiques --
655	!======================================================================
656
657	DO it=1,nbtr
658	DO k=1,klev
659	DO i=1,klon
660	d_tr_th(i,k,it)=0.
661	tr_seri(i,k,it)=MAX(tr_seri(i,k,it),0.)
662	tr_seri(i,k,it)=MIN(tr_seri(i,k,it),1.e10)
663	END DO
664	END DO
665	END DO
666
667	IF (iflag_thermals.GT.0.AND.iflag_the_trac>0) THEN
668
669	DO it=1, nbtr
670
671	CALL thermcell_dq(klon,klev,1,pdtphys,fm_therm,entr_therm, &
672	zmasse,tr_seri(1:klon,1:klev,it), &
673	d_tr_th(1:klon,1:klev,it),ztra_th,0 )
674
675	DO k=1,klev
676	DO i=1,klon
677	d_tr_th(i,k,it)=pdtphys*d_tr_th(i,k,it)
678	tr_seri(i,k,it)=MAX(tr_seri(i,k,it)+d_tr_th(i,k,it),0.)
679	END DO
680	END DO
681
682	END DO ! it
683
684	END IF ! Thermiques
685
686	!======================================================================
687	! -- Calcul de l'effet de la couche limite --
688	!======================================================================
689
690	IF (iflag_vdf_trac==1) THEN
691
692	! Injection during BL mixing
693	!
694	DO it=1, nbtr
695	!
696	IF( pbl_flg(it) /= 0 ) THEN
697	!
698	CALL cltrac(pdtphys, coefh,t_seri, &
699	tr_seri(:,:,it), source(:,it), &
700	paprs, pplay, delp, &
701	d_tr_cl(:,:,it),d_tr_dry(:,it),flux_tr_dry(:,it))
702	!
703	tr_seri(:,:,it)=tr_seri(:,:,it)+d_tr_cl(:,:,it)
704	!
705	END IF
706	!
707	END DO
708	!
709	ELSE IF (iflag_vdf_trac==0) THEN
710	!
711	! Injection of source in the first model layer
712	!
713	DO it=1,nbtr
714	d_tr_cl(:,1,it)=source(:,it)RG/delp(:,1)pdtphys
715	tr_seri(:,1,it)=tr_seri(:,1,it)+d_tr_cl(:,1,it)
716	ENDDO
717	d_tr_cl(:,2:klev,1:nbtr)=0.
718	!
719	ELSE IF (iflag_vdf_trac==-1) THEN
720	!
721	! Nothing happens
722	!
723	d_tr_cl=0.
724	!
725	ELSE
726	!
727	CALL abort_gcm('iflag_vdf_trac', 'cas non prevu',1)
728	!
729	END IF ! couche limite
730
731	!======================================================================
732	! Calcul de l'effet de la precipitation grande echelle
733	! POUR INCA le lessivage est fait directement dans INCA
734	!======================================================================
735
736	IF (lessivage) THEN
737
738	ql_incloud_ref = 10.e-4
739	ql_incloud_ref = 5.e-4
740
741
742	! calcul du contenu en eau liquide au sein du nuage
743	ql_incl = ql_incloud_ref
744	! choix du lessivage
745	!
746	IF (iflag_lscav .EQ. 3 .OR. iflag_lscav .EQ. 4) THEN
747	! ******** Olivier Boucher version (3) possibly with modified ql_incl (4)
748	!
749	DO it = 1, nbtr
750	! incloud scavenging and removal by large scale rain ! orig : ql_incl was replaced by 0.5e-3 kg/kg
751	! the value 0.5e-3 kg/kg is from Giorgi and Chameides (1986), JGR
752	! Liu (2001) proposed to use 1.5e-3 kg/kg
753
754	CALL lsc_scav(pdtphys,it,iflag_lscav,ql_incl,prfl,psfl,rneb,beta_fisrt, &
755	beta_v1,pplay,paprs,t_seri,tr_seri,d_tr_insc, &
756	d_tr_bcscav,d_tr_evapls,qPrls)
757
758	!large scale scavenging tendency
759	DO k = 1, klev
760	DO i = 1, klon
761	d_tr_ls(i,k,it)=d_tr_insc(i,k,it)+d_tr_bcscav(i,k,it)+d_tr_evapls(i,k,it)
762	tr_seri(i,k,it)=tr_seri(i,k,it)+d_tr_ls(i,k,it)
763	ENDDO
764	ENDDO
765	CALL minmaxqfi(tr_seri(:,:,it),0.,1.e33,'lsc scav it = '//solsym(it))
766	END DO !tr
767
768	ELSE IF (iflag_lscav .EQ. 2) THEN ! frac_impa, frac_nucl
769	! ********* modified old version
770
771	d_tr_lessi_nucl(:,:,:) = 0.
772	d_tr_lessi_impa(:,:,:) = 0.
773	flestottr(:,:,:) = 0.
774	! Tendance des aerosols nuclees et impactes
775	DO it = 1, nbtr
776	IF (aerosol(it)) THEN
777	his_dh(:)=0.
778	DO k = 1, klev
779	DO i = 1, klon
780	!PhH
781	zrho(i,k)=pplay(i,k)/t_seri(i,k)/RD
782	zdz(i,k)=(paprs(i,k)-paprs(i,k+1))/zrho(i,k)/RG
783	!
784	END DO
785	END DO
786
787	DO k=klev-1, 1, -1
788	DO i=1, klon
789	! d_tr_ls(i,k,it)=tr_seri(i,k,it)(frac_impa(i,k)frac_nucl(i,k)-1.)
790	dx=d_tr_ls(i,k,it)
791	his_dh(i)=his_dh(i)-dxzrho(i,k)zdz(i,k)/pdtphys ! kg/m2/s
792	evaplsc = prfl(i,k) - prfl(i,k+1) + psfl(i,k) - psfl(i,k+1)
793	! Evaporation Partielle -> Liberation Partielle 0.5*evap
794	IF ( evaplsc .LT.0..and.abs(prfl(i,k+1)+psfl(i,k+1)).gt.1.e-10) THEN
795	evaplsc = (-evaplsc)/(prfl(i,k+1)+psfl(i,k+1))
796	! evaplsc est donc positif, his_dh(i) est positif
797	!--------------
798	d_tr_evapls(i,k,it)=0.5evaplsc(d_tr_lessi_nucl(i,k+1,it) &
799	+d_tr_lessi_impa(i,k+1,it))
800	!------------- d_tr_evapls(i,k,it)=-0.5evaplsc(d_tr_lsc(i,k+1,it))
801	beta=0.5*evaplsc
802	if ((prfl(i,k)+psfl(i,k)).lt.1.e-10) THEN
803	beta=1.0*evaplsc
804	endif
805	dx=betahis_dh(i)/zrho(i,k)/zdz(i,k)pdtphys
806	his_dh(i)=(1.-beta)*his_dh(i) ! tracer from
807	d_tr_evapls(i,k,it)=dx
808	ENDIF
809	d_tr_ls(i,k,it)=tr_seri(i,k,it)(frac_impa(i,k)frac_nucl(i,k)-1.) &
810	+d_tr_evapls(i,k,it)
811
812	!--------------
813	d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) + &
814	( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
815	d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) + &
816	( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
817	!
818	! Flux lessivage total
819	flestottr(i,k,it) = flestottr(i,k,it) - &
820	( d_tr_lessi_nucl(i,k,it) + &
821	d_tr_lessi_impa(i,k,it) ) * &
822	( paprs(i,k)-paprs(i,k+1) ) / &
823	(RG * pdtphys)
824	!! Mise a jour des traceurs due a l'impaction,nucleation
825	! tr_seri(i,k,it)=tr_seri(i,k,it)frac_impa(i,k)frac_nucl(i,k)
826	!! calcul de la tendance liee au lessivage stratiforme
827	! d_tr_ls(i,k,it)=tr_seri(i,k,it)*&
828	! (1.-1./(frac_impa(i,k)*frac_nucl(i,k)))
829	!--------------
830	END DO
831	END DO
832	END IF
833	END DO
834	! ********* end modified old version
835
836	ELSE IF (iflag_lscav .EQ. 1) THEN ! frac_impa, frac_nucl
837	! ********* old version
838
839	d_tr_lessi_nucl(:,:,:) = 0.
840	d_tr_lessi_impa(:,:,:) = 0.
841	flestottr(:,:,:) = 0.
842	!=========================
843	! LESSIVAGE LARGE SCALE :
844	!=========================
845
846	! Tendance des aerosols nuclees et impactes
847	! -----------------------------------------
848	DO it = 1, nbtr
849	IF (aerosol(it)) THEN
850	DO k = 1, klev
851	DO i = 1, klon
852	d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) + &
853	( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
854	d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) + &
855	( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
856
857	!
858	! Flux lessivage total
859	! ------------------------------------------------------------
860	flestottr(i,k,it) = flestottr(i,k,it) - &
861	( d_tr_lessi_nucl(i,k,it) + &
862	d_tr_lessi_impa(i,k,it) ) * &
863	( paprs(i,k)-paprs(i,k+1) ) / &
864	(RG * pdtphys)
865	!
866	! Mise a jour des traceurs due a l'impaction,nucleation
867	! ----------------------------------------------------------------------
868	tr_seri(i,k,it)=tr_seri(i,k,it)frac_impa(i,k)frac_nucl(i,k)
869	END DO
870	END DO
871	END IF
872	END DO
873
874	! ********* end old version
875	ENDIF ! iflag_lscav . EQ. 1, 2, 3 or 4
876	!
877	END IF ! lessivage
878
879	!=============================================================
880	! Ecriture des sorties
881	!=============================================================
882	#ifdef CPP_IOIPSL
883	! INCLUDE "write_histrac.h"
884	#endif
885
886	END SUBROUTINE phytrac
887
888	END MODULE

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

Download in other formats:

Original Format