Context Navigation

phytrac.F @ 1134

Last change on this file since 1134 was 1114, checked in by jghattas, 16 years ago

Creation du module infotrac:

contient les variables de advtrac.h
contient la subroutine iniadvtrac renommer en infotrac_init
le nombre des traceurs est lu dans tracer.def en dynamique (ou par default ou recu par INCA)
ce module est utilise dans la dynamique et la physique
contient aussi la variable nbtr qui avant etait stockee dans dimphy

Le fichier advtrac.h n'existe plus.
La compilation ne prend plus en compte le nombre de traceur.

/JG

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 26.5 KB

Line
1	!
2	c
3	c
4	SUBROUTINE phytrac (rnpb,
5	I nstep,
6	I julien,
7	I gmtime,
8	I debutphy,
9	I lafin,
10	I nlon,
11	I nlev,
12	I pdtphys,
13	I u,
14	I v,
15	I t_seri,
16	I paprs,
17	I pplay,
18	I pmfu,
19	I pmfd,
20	I pen_u,
21	I pde_u,
22	I pen_d,
23	I pde_d,
24	I cdragh,
25	I coefh,
26	I fm_therm,
27	I entr_therm,
28	I yu1,
29	I yv1,
30	I ftsol,
31	I pctsrf,
32	I xlat,
33	I frac_impa,
34	I frac_nucl,
35	I xlon,
36	I presnivs,
37	I pphis,
38	I pphi,
39	I albsol,
40	I sh,
41	I rh,
42	I cldfra,
43	I rneb,
44	I diafra,
45	I cldliq,
46	I itop_con,
47	I ibas_con,
48	I pmflxr,
49	I pmflxs,
50	I prfl,
51	I psfl,
52	I da,
53	I phi,
54	I mp,
55	I upwd,
56	I dnwd,
57	I aerosol_couple,
58	I flxmass_w,
59	I tau_inca,
60	I piz_inca,
61	I cg_inca,
62	I ccm,
63	I rfname,
64	O tr_seri)
65
66	USE ioipsl
67	USE dimphy
68	USE infotrac
69	USE mod_grid_phy_lmdz
70	USE mod_phys_lmdz_para
71	USE comgeomphy
72	USE iophy
73	USE vampir
74
75	IMPLICIT none
76	c======================================================================
77	c Auteur(s) FH
78	c Objet: Moniteur general des tendances traceurs
79	c
80	cAA Remarques en vrac:
81	cAA--------------------
82	cAA 2/ Le choix du radon et du pb se fait juste avec un data
83	cAA (peu propre). Peut-etre pourrait-on prevoir dans l'avenir
84	cAA une variable "type de traceur"
85	c======================================================================
86	#include "YOMCST.h"
87	#include "dimensions.h"
88	#include "indicesol.h"
89	#include "clesphys.h"
90	#include "temps.h"
91	#include "paramet.h"
92	#include "control.h"
93	#include "thermcell.h"
94	c======================================================================
95
96	c Arguments:
97	c
98	c EN ENTREE:
99	c ==========
100	c
101	c divers:
102	c -------
103	c
104	integer nlon ! nombre de points horizontaux
105	integer nlev ! nombre de couches verticales
106	integer nstep ! appel physique
107	integer julien !jour julien
108	integer itop_con(nlon)
109	integer ibas_con(nlon)
110	real gmtime
111	real pdtphys ! pas d'integration pour la physique (seconde)
112	real t_seri(nlon,nlev) ! temperature
113	real tr_seri(nlon,nlev,nbtr) ! traceur
114	real u(klon,klev)
115	real v(klon,klev)
116	real sh(nlon,nlev) ! humidite specifique
117	real rh(nlon,nlev) ! humidite relative
118	real cldliq(nlon,nlev) ! eau liquide nuageuse
119	real cldfra(nlon,nlev) ! fraction nuageuse (tous les nuages)
120	real diafra(nlon,nlev) ! fraction nuageuse (convection ou stratus artificiels)
121	real rneb(nlon,nlev) ! fraction nuageuse (grande echelle)
122	real albsol(nlon) ! albedo surface
123	real paprs(nlon,nlev+1) ! pression pour chaque inter-couche (en Pa)
124	real ps(nlon) ! pression surface
125	real pplay(nlon,nlev) ! pression pour le mileu de chaque couche (en Pa)
126	real pphi(nlon,klev) ! geopotentiel
127	real pphis(klon)
128	REAL presnivs(klev)
129	logical debutphy ! le flag de l'initialisation de la physique
130	logical lafin ! le flag de la fin de la physique
131	c Olivia
132	integer nsplit
133	REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1) !--lessivage convection
134	REAL prfl(klon,klev+1), psfl(klon,klev+1) !--lessivage large-scale
135	LOGICAL aerosol_couple
136
137	REAL flxmass_w(klon,klev)
138	CHARACTER(len=8) :: solsym(nbtr)
139	integer la
140	REAL :: tau_inca(klon,klev,9,2)
141	REAL :: piz_inca(klon,klev,9,2)
142	REAL :: cg_inca(klon,klev,9,2)
143	character*4 :: rfname(9)
144	REAL :: ccm(klon,klev,2)
145	c
146	c convection:
147	c -----------
148	c
149	REAL pmfu(nlon,nlev) ! flux de masse dans le panache montant
150	REAL pmfd(nlon,nlev) ! flux de masse dans le panache descendant
151	REAL pen_u(nlon,nlev) ! flux entraine dans le panache montant
152
153	c
154	c thermiques:
155	c -----------
156	c
157	real fm_therm(klon,klev+1),entr_therm(klon,klev)
158	real fm_therm1(klon,klev)
159	c
160	REAL pde_u(nlon,nlev) ! flux detraine dans le panache montant
161	REAL pen_d(nlon,nlev) ! flux entraine dans le panache descendant
162	REAL pde_d(nlon,nlev) ! flux detraine dans le panache descendant
163	c KE
164	real da(nlon,nlev),phi(nlon,nlev,nlev),mp(nlon,nlev)
165	REAL upwd(nlon,nlev) ! saturated updraft mass flux
166	REAL dnwd(nlon,nlev) ! saturated downdraft mass flux
167
168	c
169	c Couche limite:
170	c --------------
171	c
172	REAL cdragh(nlon,nlev)! coeff drag pour T et Q
173	REAL coefh(nlon,nlev) ! coeff melange CL
174	REAL yu1(nlon) ! vents au premier niveau
175	REAL yv1(nlon) ! vents au premier niveau
176	REAL xlat(nlon) ! latitudes pour chaque point
177	REAL xlon(nlon) ! longitudes pour chaque point
178
179	c
180	c Lessivage:
181	c ----------
182	c
183	c pour le ON-LINE
184	c
185	REAL frac_impa(nlon,nlev) ! fraction d'aerosols impactes
186	REAL frac_nucl(nlon,nlev) ! fraction d'aerosols nuclees
187	c
188	cAA
189	cAA Arguments necessaires pour les sources et puits de traceur:
190	cAA ----------------
191	cAA
192	real ftsol(nlon,nbsrf) ! Temperature du sol (surf)(Kelvin)
193	real pctsrf(nlon,nbsrf) ! Pourcentage de sol f(nature du sol)
194	c abder
195	real pftsol1(nlon),pftsol2(nlon),pftsol3(nlon),pftsol4(nlon)
196	real ppsrf1(nlon),ppsrf2(nlon),ppsrf3(nlon),ppsrf4(nlon)
197	c fin
198	cAA ----------------------------
199	cAA VARIABLES LOCALES TRACEURS
200	cAA ----------------------------
201	cAA
202	cAA Sources et puits des traceurs:
203	cAA ------------------------------
204	cAA
205	cAA Pour l'instant seuls les cas du rn et du pb ont ete envisages.
206
207	REAL source(klon,nbtr) ! a voir lorsque le flux est prescrit
208	cAA
209	cAA Pour la source de radon et son reservoir de sol
210	cAA ................................................
211
212	REAL,save,allocatable :: trs(:,:) ! Conc. radon ds le sol
213	c$OMP THREADPRIVATE(trs)
214	REAL :: trs_tmp(klon_glo)
215
216	REAL,save,allocatable :: masktr(:,:) ! Masque reservoir de sol traceur
217	c Masque de l'echange avec la surface
218	c (1 = reservoir) ou (possible => 1 )
219	c$OMP THREADPRIVATE(masktr)
220	REAL,save,allocatable :: fshtr(:,:) ! Flux surfacique dans le reservoir de sol
221	c$OMP THREADPRIVATE(fshtr)
222	REAL,save,allocatable :: hsoltr(:) ! Epaisseur equivalente du reservoir de sol
223	c$OMP THREADPRIVATE(hsoltr)
224	REAL,save,allocatable :: tautr(:) ! Constante de decroissance radioactive
225	c$OMP THREADPRIVATE(tautr)
226	REAL,save,allocatable :: vdeptr(:) ! Vitesse de depot sec dans la couche Brownienne
227	c$OMP THREADPRIVATE(vdeptr)
228	REAL,save,allocatable :: scavtr(:) ! Coefficient de lessivage
229	c$OMP THREADPRIVATE(scavtr)
230	cAA
231	CHARACTER*2 itn
232	C maf ioipsl
233	CHARACTER*2 str2
234	INTEGER nhori, nvert
235	REAL zsto, zout, zjulian
236	INTEGER, SAVE :: nid_tra
237	c$OMP THREADPRIVATE(nid_tra)
238	INTEGER, SAVE :: nid_tra2,nid_tra3
239	c$OMP THREADPRIVATE(nid_tra2,nid_tra3)
240	INTEGER ndex(1)
241	INTEGER ndex2d(iim(jjm+1)),ndex3d(iim(jjm+1)*klev)
242	REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
243	REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
244	REAL zx_tmp_2d(iim,jjm+1), zx_tmp_3d(iim,jjm+1,klev)
245	REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1)
246	c
247	integer itau_w ! pas de temps ecriture = nstep + itau_phy
248	c
249	logical ok_sync
250	parameter (ok_sync = .true.)
251	C
252	C nature du traceur
253	c
254	logical,save,allocatable :: aerosol(:) ! Nature du traceur
255	c ! aerosol(it) = true => aerosol
256	c ! aerosol(it) = false => gaz
257	c ! nat_trac(it) = 1. aerosol
258	logical,save,allocatable :: clsol(:) ! clsol(it) = true => CL sol calculee
259	logical,save,allocatable :: radio(:) ! radio(it)=true => decroisssance radioactive
260	c$OMP THREADPRIVATE(aerosol,clsol,radio)
261	C
262	c======================================================================
263	c
264	c Declaration des procedures appelees
265	c
266	c--modif convection tiedtke
267	INTEGER i, k, it
268	INTEGER iq, iiq
269	REAL delp(klon,klev)
270	c--end modif
271	c
272	c Variables liees a l'ecriture de la bande histoire physique
273	c
274	c Variables locales pour effectuer les appels en serie
275	c----------------------------------------------------
276	c
277	REAL d_tr(klon,klev), d_trs(klon) ! tendances de traceurs
278	REAL d_tr_cl(klon,klev,nbtr) ! tendance de traceurs couche limite
279	REAL d_tr_cv(klon,klev,nbtr) ! tendance de traceurs conv pour chq traceur
280	REAL d_tr_th(klon,klev,nbtr) ! la tendance des thermiques
281	REAL d_tr_dec(klon,klev,nbtr) ! la tendance de la decroissance
282	c ! radioactive du rn - > pb
283	REAL d_tr_lessi_impa(klon,klev,nbtr) ! la tendance du lessivage
284	c ! par impaction
285	REAL d_tr_lessi_nucl(klon,klev,nbtr) ! la tendance du lessivage
286	c ! par nucleation
287	REAL fluxrn(klon,klev)
288	REAL fluxpb(klon,klev)
289	REAL pbimpa(klon,klev)
290	REAL pbnucl(klon,klev)
291	REAL rn(klon,klev)
292	REAL pb(klon,klev)
293	REAL flestottr(klon,klev,nbtr) ! flux de lessivage
294	c ! dans chaque couche
295	real zmasse(klon,klev)
296	real ztra_th(klon,klev)
297
298	C
299	character*20 modname
300	character*80 abort_message
301	c
302	c Controles
303	c-------------
304	logical first,couchelimite,convection,lessivage,sorties,
305	s rnpb,inirnpb
306	save first,couchelimite,convection,lessivage,
307	s sorties,inirnpb
308	c$OMP THREADPRIVATE(first,couchelimite,convection,lessivage,
309	c$OMP+ sorties,inirnpb)
310	c data first,couchelimite,convection,lessivage,sorties
311	c s /.true.,.true.,.false.,.true.,.true./
312	c Olivia
313	data first,couchelimite,convection,lessivage,
314	s sorties
315	s /.true.,.true.,.true.,.true.,.true./
316
317	! Variables needed for configuration with INCA
318	INTEGER :: lastgas
319	INTEGER :: ncsec
320
321	REAL, PARAMETER :: dry_mass = 28.966
322	REAL, POINTER :: hbuf(:)
323	REAL, ALLOCATABLE :: obuf(:)
324	REAL :: calday
325	REAL :: pdel(klon,klev)
326	c
327	c======================================================================
328
329	modname='phytrac'
330
331	ps(:)=paprs(:,1)
332
333	if (debutphy) then
334	allocate( trs(klon,nbtr) )
335	allocate( masktr(klon,nbtr))
336	allocate( fshtr(klon,nbtr) )
337	allocate( hsoltr(nbtr))
338	allocate( tautr(nbtr))
339	allocate( vdeptr(nbtr))
340	allocate( scavtr(nbtr))
341	allocate( aerosol(nbtr))
342	allocate( clsol(nbtr))
343	allocate( radio(nbtr))
344
345
346	! FH 2008/05/09 correction de la frequence d'ecriture des traceurs
347	! ecrit_tra = FLOAT(NINT(86400./pdtphys *ecritphy))
348	print*,'dans phytrac ',pdtphys,ecrit_tra
349
350	inirnpb=rnpb
351	PRINT*, 'La frequence de sortie traceurs est ', ecrit_tra
352	C
353	c=============================================================
354	c Initialisation des sorties
355	c=============================================================
356
357	#ifdef CPP_IOIPSL
358	#include "ini_histrac.h"
359	#endif
360
361	c======================================================================
362	c Initialisation de certaines variables pour le Rn et le Pb
363	c======================================================================
364
365	c Initialisation du traceur dans le sol (couche limite radonique)
366	c
367	c print*,'valeur de debut dans phytrac :',debutphy
368	trs(:,:) = 0.
369	c$OMP MASTER
370	if (is_mpi_root) then
371	trs_tmp(:)=0.
372	open (99,file='starttrac',status='old',
373	. err=999,form='formatted')
374	read(99,*) (trs_tmp(i),i=1,klon_glo)
375	999 close(99)
376	endif
377	c$OMP END MASTER
378	call Scatter(trs_tmp,trs(:,1))
379
380	c print*, 'apres starttrac'
381
382	c Initialisation de la fraction d'aerosols lessivee
383	c
384	d_tr_lessi_impa(:,:,:) = 0.
385	d_tr_lessi_nucl(:,:,:) = 0.
386	c
387	c Initialisation de la nature des traceurs
388	c
389	DO it = 1, nbtr
390	aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut
391	radio(it) = .FALSE. ! Par defaut pas de passage par radiornpb
392	clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit
393	ENDDO
394	c
395	ENDIF ! fin debutphy
396	c Initialisation du traceur dans le sol (couche limite radonique)
397	if(inirnpb) THEN
398	c
399	radio(1)= .true.
400	radio(2)= .true.
401	clsol(1)= .true.
402	clsol(2)= .true.
403	aerosol(2) = .TRUE. ! le Pb est un aerosol
404	c
405	call initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr
406	. ,vdeptr,scavtr)
407	inirnpb=.false.
408	endif
409
410
411	IF (config_inca == 'none') THEN
412	DO i=1,nlon
413	pftsol1(i) = ftsol(i,1)
414	pftsol2(i) = ftsol(i,2)
415	pftsol3(i) = ftsol(i,3)
416	pftsol4(i) = ftsol(i,4)
417
418	ppsrf1(i) = pctsrf(i,1)
419	ppsrf2(i) = pctsrf(i,2)
420	ppsrf3(i) = pctsrf(i,3)
421	ppsrf4(i) = pctsrf(i,4)
422
423	ENDDO
424
425	ELSE ! config_inca /=none
426	#ifdef INCA
427	call VTe(VTphysiq)
428	call VTb(VTinca)
429
430	!======================================================================
431	! Chimie
432	!======================================================================
433
434	calday = FLOAT(julien) + gmtime
435	ncsec = NINT (86400.*gmtime)
436
437	DO k = 1, nlev
438	pdel(:,k) = paprs(:,k) - paprs (:,k+1)
439	END DO
440
441	IF (config_inca == 'aero') THEN
442	CALL aerosolmain (aerosol_couple,
443	$ tr_seri,
444	$ pdtphys,
445	$ pplay,
446	$ pdel,
447	$ prfl,
448	$ pmflxr,
449	$ psfl,
450	$ pmflxs,
451	$ pmfu,
452	$ itop_con,
453	$ ibas_con,
454	$ pphi,
455	$ airephy, ! paire,
456	$ nstep,
457	$ rneb, ! for chimiaq
458	$ t_seri, ! for chimiaq
459	$ rh,
460	$ tau_inca,
461	$ piz_inca,
462	$ cg_inca,
463	$ rfname,
464	$ ccm,
465	$ lafin)
466	END IF
467
468	CALL chemmain (tr_seri, !mmr
469	$ nstep, !nstep
470	$ calday, !calday
471	$ julien, !ncdate
472	$ ncsec, !ncsec
473	$ 1, !lat
474	$ pdtphys, !delt
475	$ paprs(1,1), !ps
476	$ pplay, !pmid
477	$ pdel, !pdel
478	$ airephy,
479	$ pctsrf(1,1),!oro
480	$ ftsol, !tsurf
481	$ albsol, !albs
482	$ pphi, !zma
483	$ pphis, !phis
484	$ cldfra, !cldfr
485	$ rneb, !cldfr_st
486	$ diafra, !cldfr_cv
487	$ itop_con, !cldtop
488	$ ibas_con, !cldbot
489	$ cldliq, !cwat
490	$ prfl, !flxrst
491	$ pmflxr, !flxrcv
492	$ psfl, !flxsst
493	$ pmflxs, !flxscv
494	$ pmfu, !flxupd
495	$ flxmass_w, !flxmass_w
496	$ t_seri, !tfld
497	$ sh, !sh
498	$ rh, !rh
499	$ .false., !wrhstts
500	$ hbuf, !hbuf
501	$ obuf, !obuf
502	$ iip1, !nx
503	$ jjp1, !ny
504	$ source,
505	$ solsym)
506
507
508	call VTe(VTinca)
509	call VTb(VTphysiq)
510	#endif
511	END IF ! config_inca
512	c======================================================================
513	c Calcul de l'effet de la convection
514	c======================================================================
515	c print*,'Avant convection'
516	do it=1,nbtr
517	WRITE(itn,'(i2)') it
518	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
519	enddo
520
521	if (convection) then
522
523	c print*,'Pas de temps dans phytrac : ',pdtphys
524	DO it=1, nbtr
525
526	IF ( config_inca/='none') THEN
527	IF ( conv_flg(it) == 0 ) CYCLE
528	END IF
529
530	if (iflag_con.lt.2) then
531	d_tr_cv=0.
532	else if (iflag_con.eq.2) then
533	c tiedke
534	CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
535	. pplay, paprs, tr_seri(1,1,it), d_tr_cv(1,1,it))
536	else
537	c KE
538	call cvltr(pdtphys, da, phi, mp, paprs,pplay, tr_seri(1,1,it),
539	. upwd,dnwd,d_tr_cv(1,1,it))
540	endif
541
542	DO k = 1, nlev
543	DO i = 1, klon
544	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it)
545	ENDDO
546	ENDDO
547
548	IF (config_inca == 'none') THEN
549	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '//itn)
550	ELSE
551	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '
552	. //solsym(it))
553	END IF
554
555	ENDDO
556
557	endif ! convection
558	c print*,'Apres convection'
559	c do it=1,nbtr
560	c WRITE(itn,'(i1)') it
561	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
562	c enddo
563
564
565	c======================================================================
566	c Calcul de l'effet des thermiques
567	c======================================================================
568
569	do k=1,klev
570	do i=1,klon
571	zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
572	enddo
573	enddo
574
575	c print*,'masse dans ph ',zmasse
576	do it=1,nbtr
577	do k=1,klev
578	do i=1,klon
579	d_tr_th(i,k,it)=0.
580	tr_seri(i,k,it)=max(tr_seri(i,k,it),0.)
581	tr_seri(i,k,it)=min(tr_seri(i,k,it),1.e10)
582	enddo
583	enddo
584	enddo
585
586	if (iflag_thermals.gt.0) then
587	c print*,'calcul de leffet des thermiques'
588	nsplit=10
589	DO it=1, nbtr
590	c WRITE(itn,'(i1)') it
591	c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'conv it='//itn)
592	c print*,'avant dqthermiquesretro'
593	c call dump2d(iim,jjm-1,tr_seri(2,1,1),'TR_SERI ')
594
595	do isplit=1,nsplit
596	c Abderr 25 11 02
597	C Thermiques
598	c print*,'Avant dans phytrac'
599	call dqthermcell(klon,klev,pdtphys/nsplit
600	. ,fm_therm,entr_therm,zmasse
601	. ,tr_seri(1:klon,1:klev,it),d_tr,ztra_th)
602
603	do k=1,klev
604	do i=1,klon
605	d_tr(i,k)=pdtphys*d_tr(i,k)/nsplit
606	d_tr_th(i,k,it)=d_tr_th(i,k,it)+d_tr(i,k)
607	tr_seri(i,k,it)=max(tr_seri(i,k,it)+d_tr(i,k),0.)
608	enddo
609	enddo
610	enddo ! nsplit
611	c print*,'apres thermiques'
612	c call dump2d(iim,jjm-1,d_tr_th(1,1,1),'d_tr_th ')
613	c do k=1,klev
614	c print*,'d_tr_th(',k,')=',tr_seri(280,k,1)
615	c enddo
616
617	c WRITE(itn,'(i1)') it
618	c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'therm it='//itn)
619	ENDDO ! it
620	endif ! Thermiques
621	c print*,'ATTENTION: sdans thermniques'
622
623	c======================================================================
624	c Calcul de l'effet de la couche limite
625	c======================================================================
626	c print *,'Avant couchelimite'
627	c do it=1,nbtr
628	c WRITE(itn,'(i1)') it
629	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
630	c enddo
631
632	if (couchelimite) then
633
634	DO k = 1, nlev
635	DO i = 1, klon
636	delp(i,k) = paprs(i,k)-paprs(i,k+1)
637	ENDDO
638	ENDDO
639
640	C maf modif pour tenir compte du cas rnpb + traceur
641	DO it=1, nbtr
642
643	IF ( config_inca/='none' ) THEN
644	IF( pbl_flg(it) == 0 ) CYCLE
645	END IF
646
647	c print *,'it',it,clsol(it)
648	if (clsol(it)) then ! couche limite avec quantite dans le sol calculee
649	CALL cltracrn(it, pdtphys, yu1, yv1,
650	e cdragh, coefh,t_seri,ftsol,pctsrf,
651	e tr_seri(1,1,it),trs(1,it),
652	e paprs, pplay, delp,
653	e masktr(1,it),fshtr(1,it),hsoltr(it),
654	e tautr(it),vdeptr(it),
655	e xlat,
656	s d_tr_cl(1,1,it),d_trs)
657	DO k = 1, nlev
658	DO i = 1, klon
659	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
660	ENDDO
661	ENDDO
662	c
663	c Traceur ds sol
664	c
665	DO i = 1, klon
666	trs(i,it) = trs(i,it) + d_trs(i)
667	END DO
668	C
669	C maf provisoire suppression des prints
670	C WRITE(itn,'(i1)') it
671	C CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
672	else ! couche limite avec flux prescrit
673	#ifndef INCA
674
675	Cmaf provisoire source / traceur a creer
676	DO i=1, klon
677	source(i,it) = 0.0 ! pas de source, pour l'instant
678	ENDDO
679	C
680	#endif
681	CALL cltrac(pdtphys, coefh,t_seri,
682	s tr_seri(1,1,it), source(:,it),
683	e paprs, pplay, delp,
684	s d_tr_cl(1,1,it))
685	DO k = 1, nlev
686	DO i = 1, klon
687	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
688	ENDDO
689	ENDDO
690	Cmaf WRITE(itn,'(i1)') it
691	cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
692	endif
693	ENDDO
694	c
695	endif ! couche limite
696
697	c print*,'Apres couchelimite'
698	c do it=1,nbtr
699	c WRITE(itn,'(i1)') it
700	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
701	c enddo
702
703	c======================================================================
704	c Calcul de l'effet du puits radioactif
705	c======================================================================
706
707	C MAF il faudrait faire une modification pour passer dans radiornpb
708	c si radio=true mais pour l'instant radiornpb propre au cas rnpb
709	if(rnpb) then
710	c print *, 'decroissance radiactive activee'
711	call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
712	C
713	DO it=1,nbtr
714	if(radio(it)) then
715	DO k = 1, nlev
716	DO i = 1, klon
717	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
718	ENDDO
719	ENDDO
720	WRITE(itn,'(i1)') it
721	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn)
722	endif
723	ENDDO
724	c
725	endif ! rnpb decroissance radioactive
726	C
727	c======================================================================
728	c Calcul de l'effet de la precipitation
729	c======================================================================
730
731	c print*,'LESSIVAGE =',lessivage
732	IF (lessivage) THEN
733
734	c print*,'avant lessivage'
735
736	d_tr_lessi_nucl(:,:,:) = 0.
737	d_tr_lessi_impa(:,:,:) = 0.
738	flestottr(:,:,:) = 0.
739	c
740	c tendance des aerosols nuclees et impactes
741	c
742	DO it = 1, nbtr
743	IF (aerosol(it)) THEN
744	DO k = 1, nlev
745	DO i = 1, klon
746	d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +
747	s ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
748	d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +
749	s ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
750	ENDDO
751	ENDDO
752	ENDIF
753	ENDDO
754	c
755	c Mises a jour des traceurs + calcul des flux de lessivage
756	c Mise a jour due a l'impaction et a la nucleation
757	c
758	c call dump2d(iim,jjm-1,frac_impa(2:klon-1,10),'FRACIMPA')
759	c call dump2d(iim,jjm-1,frac_nucl(2:klon-1,10),'FRACNUCL')
760	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3')
761	DO it = 1, nbtr
762	c print*,'IT=',it,aerosol(it)
763	IF (aerosol(it)) THEN
764	c print*,'IT=',it,' On lessive'
765	DO k = 1, nlev
766	DO i = 1, klon
767	tr_seri(i,k,it)=tr_seri(i,k,it)
768	s frac_impa(i,k)frac_nucl(i,k)
769	ENDDO
770	ENDDO
771	ENDIF
772	ENDDO
773	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3B')
774	c
775	c Flux lessivage total
776	c
777	DO it = 1, nbtr
778	DO k = 1, nlev
779	DO i = 1, klon
780	flestottr(i,k,it) = flestottr(i,k,it) -
781	s ( d_tr_lessi_nucl(i,k,it) +
782	s d_tr_lessi_impa(i,k,it) ) *
783	s ( paprs(i,k)-paprs(i,k+1) ) /
784	s (RG * pdtphys)
785	ENDDO
786	ENDDO
787	c
788	Cmaf WRITE(itn,'(i1)') it
789	Cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
790	ENDDO
791	c
792	c print*,'apres lessivage'
793	ENDIF
794	Cc
795	DO k = 1, nlev
796	DO i = 1, klon
797	fluxrn(i,k) = flestottr(i,k,1)
798	fluxpb(i,k) = flestottr(i,k,2)
799	rn(i,k) = tr_seri(i,k,1)
800	pb(i,k) = tr_seri(i,k,2)
801	pbnucl(i,k)=d_tr_lessi_nucl(i,k,2)
802	pbimpa(i,k)=d_tr_lessi_impa(i,k,2)
803	ENDDO
804	ENDDO
805
806	c=============================================================
807	c Ecriture des sorties
808	c=============================================================
809
810	#ifdef CPP_IOIPSL
811	#include "write_histrac.h"
812	#endif
813
814	c=============================================================
815
816	if (lafin) then
817	print*, 'c est la fin de la physique'
818	call Gather(trs(:,1),trs_tmp)
819	c$OMP MASTER
820	if (is_mpi_root) then
821
822	open (99,file='restarttrac', form='formatted')
823	do i=1,klon_glo
824	write(99,*) trs_tmp(i)
825	enddo
826	PRINT*, 'Ecriture du fichier restarttrac'
827	close(99)
828	endif
829	c$OMP END MASTER
830	else
831	c print*, 'physique pas fini'
832	endif
833
834
835	RETURN
836	END

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

Download in other formats:

Original Format