Context Navigation

phytrac.F @ 970

Last change on this file since 970 was 959, checked in by lsce, 17 years ago
Ajoute du parametre config_inca dans conf_gcm.F config_inca='none'(sans INCA, par defaut) config_inca='chem'(avec INCA config chemie) config_inca='aero'(avec INCA config aerosol) Menage parmis les cles CPP INCA Enleve le calcul d'omega dans calfis.F et active le calcul correspondant dans advtrac.F(avant uniquement pour INCA).
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 27.0 KB

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

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

Download in other formats:

Original Format