Context Navigation

phytrac.F @ 557

Last change on this file since 557 was 541, checked in by lmdzadmin, 20 years ago

Convergence avec la version d'Olivia Coindreau incluant:

le offline
les thermiques
mellor & yamada dans la couche limite

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

File size: 25.9 KB

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