Context Navigation

phytrac_mod.F90 @ 1820

Last change on this file since 1820 was 1820, checked in by Ehouarn Millour, 11 years ago
Correction d'un bug dans phytrac affectant la compilation en openmp. Attention à l'heure actuelle il n'y a plus convergence entre la version séquentielle et la version omp. UG ..................... Bug correction in phytrac, no compiling in openmp as well as in sequential. Be warned that convergence have been lost between sequential and openmp versions. UG
File size: 28.7 KB

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

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

Download in other formats:

Original Format