Context Navigation

phytrac.F @ 1203

Last change on this file since 1203 was 1146, checked in by Laurent Fairhead, 16 years ago
Réintegration dans le tronc des modifications issues de la branche LMDZ-dev comprises entre la révision 1074 et 1145 Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes résultats entre la trunk et la dev LF
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
674	IF (config_inca == 'none') THEN
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	END IF
680
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