Context Navigation

phytrac.F @ 1105

Last change on this file since 1105 was 1067, checked in by Laurent Fairhead, 16 years ago

Modifications lie au premier niveau du modele pour la diffusion turbulent

du vent.

Preparation pour un couplage des courrant oceaniques.

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

File size: 27.1 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 nqmax,
11	I nlon,
12	I nlev,
13	I pdtphys,
14	I u,
15	I v,
16	I t_seri,
17	I paprs,
18	I pplay,
19	I pmfu,
20	I pmfd,
21	I pen_u,
22	I pde_u,
23	I pen_d,
24	I pde_d,
25	I cdragh,
26	I coefh,
27	I fm_therm,
28	I entr_therm,
29	I yu1,
30	I yv1,
31	I ftsol,
32	I pctsrf,
33	I xlat,
34	I frac_impa,
35	I frac_nucl,
36	I xlon,
37	I presnivs,
38	I pphis,
39	I pphi,
40	I albsol,
41	I sh,
42	I rh,
43	I cldfra,
44	I rneb,
45	I diafra,
46	I cldliq,
47	I itop_con,
48	I ibas_con,
49	I pmflxr,
50	I pmflxs,
51	I prfl,
52	I psfl,
53	I da,
54	I phi,
55	I mp,
56	I upwd,
57	I dnwd,
58	I aerosol_couple,
59	I flxmass_w,
60	I tau_inca,
61	I piz_inca,
62	I cg_inca,
63	I ccm,
64	I rfname,
65	O tr_seri)
66
67	USE ioipsl
68	USE dimphy
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 1/ le call phytrac se fait avec nqmax-2 donc nous avons bien
83	cAA les vrais traceurs (nbtr) dans phytrac (pas la vapeur ni eau liquide)
84	cAA 2/ Le choix du radon et du pb se fait juste avec un data
85	cAA (peu propre). Peut-etre pourrait-on prevoir dans l'avenir
86	cAA une variable "type de traceur"
87	c======================================================================
88	#include "YOMCST.h"
89	#include "dimensions.h"
90	#include "indicesol.h"
91	#include "clesphys.h"
92	#include "temps.h"
93	#include "paramet.h"
94	#include "control.h"
95	#include "advtrac.h"
96	#include "thermcell.h"
97	c======================================================================
98
99	c Arguments:
100	c
101	c EN ENTREE:
102	c ==========
103	c
104	c divers:
105	c -------
106	c
107	integer nlon ! nombre de points horizontaux
108	integer nlev ! nombre de couches verticales
109	integer nqmax ! nombre de traceurs auxquels on applique la physique
110	integer nstep ! appel physique
111	integer julien !jour julien
112	integer itop_con(nlon)
113	integer ibas_con(nlon)
114	real gmtime
115	real pdtphys ! pas d'integration pour la physique (seconde)
116	real t_seri(nlon,nlev) ! temperature
117	real tr_seri(nlon,nlev,nbtr) ! traceur
118	real u(klon,klev)
119	real v(klon,klev)
120	real sh(nlon,nlev) ! humidite specifique
121	real rh(nlon,nlev) ! humidite relative
122	real cldliq(nlon,nlev) ! eau liquide nuageuse
123	real cldfra(nlon,nlev) ! fraction nuageuse (tous les nuages)
124	real diafra(nlon,nlev) ! fraction nuageuse (convection ou stratus artificiels)
125	real rneb(nlon,nlev) ! fraction nuageuse (grande echelle)
126	real albsol(nlon) ! albedo surface
127	real paprs(nlon,nlev+1) ! pression pour chaque inter-couche (en Pa)
128	real ps(nlon) ! pression surface
129	real pplay(nlon,nlev) ! pression pour le mileu de chaque couche (en Pa)
130	real pphi(nlon,klev) ! geopotentiel
131	real pphis(klon)
132	REAL presnivs(klev)
133	logical debutphy ! le flag de l'initialisation de la physique
134	logical lafin ! le flag de la fin de la physique
135	c Olivia
136	integer nsplit
137	REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1) !--lessivage convection
138	REAL prfl(klon,klev+1), psfl(klon,klev+1) !--lessivage large-scale
139	LOGICAL aerosol_couple
140
141	REAL flxmass_w(klon,klev)
142	CHARACTER(len=8) :: solsym(nqmax)
143	integer la
144	REAL :: tau_inca(klon,klev,9,2)
145	REAL :: piz_inca(klon,klev,9,2)
146	REAL :: cg_inca(klon,klev,9,2)
147	character*4 :: rfname(9)
148	REAL :: ccm(klon,klev,2)
149	c
150	c convection:
151	c -----------
152	c
153	REAL pmfu(nlon,nlev) ! flux de masse dans le panache montant
154	REAL pmfd(nlon,nlev) ! flux de masse dans le panache descendant
155	REAL pen_u(nlon,nlev) ! flux entraine dans le panache montant
156
157	c
158	c thermiques:
159	c -----------
160	c
161	real fm_therm(klon,klev+1),entr_therm(klon,klev)
162	real fm_therm1(klon,klev)
163	c
164	REAL pde_u(nlon,nlev) ! flux detraine dans le panache montant
165	REAL pen_d(nlon,nlev) ! flux entraine dans le panache descendant
166	REAL pde_d(nlon,nlev) ! flux detraine dans le panache descendant
167	c KE
168	real da(nlon,nlev),phi(nlon,nlev,nlev),mp(nlon,nlev)
169	REAL upwd(nlon,nlev) ! saturated updraft mass flux
170	REAL dnwd(nlon,nlev) ! saturated downdraft mass flux
171
172	c
173	c Couche limite:
174	c --------------
175	c
176	REAL cdragh(nlon,nlev)! coeff drag pour T et Q
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 cdragh, 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