Context Navigation

phytrac.F @ 777

Last change on this file since 777 was 776, checked in by Laurent Fairhead, 17 years ago
Suite du merge entre la version et la HEAD: quelques modifications de Yann sur le LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 29.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	#ifdef INCA
60	I flxmass_w,
61	#endif
62	O tr_seri)
63
64	USE ioipsl
65	USE dimphy
66	USE mod_grid_phy_lmdz
67	USE mod_phys_lmdz_para
68	USE comgeomphy
69	USE iophy
70	USE vampir
71
72	IMPLICIT none
73	c======================================================================
74	c Auteur(s) FH
75	c Objet: Moniteur general des tendances traceurs
76	c
77	cAA Remarques en vrac:
78	cAA--------------------
79	cAA 1/ le call phytrac se fait avec nqmax-2 donc nous avons bien
80	cAA les vrais traceurs (nbtr) dans phytrac (pas la vapeur ni eau liquide)
81	cAA 2/ Le choix du radon et du pb se fait juste avec un data
82	cAA (peu propre). Peut-etre pourrait-on prevoir dans l'avenir
83	cAA une variable "type de traceur"
84	c======================================================================
85	#include "YOMCST.h"
86	#include "dimensions.h"
87	cym#include "dimphy.h"
88	#include "indicesol.h"
89	#include "clesphys.h"
90	#include "temps.h"
91	#include "paramet.h"
92	#include "control.h"
93	cym#include "comgeomphy.h"
94	#include "advtrac.h"
95	#include "thermcell.h"
96	c======================================================================
97
98	c Arguments:
99	c
100	c EN ENTREE:
101	c ==========
102	c
103	c divers:
104	c -------
105	c
106	integer nlon ! nombre de points horizontaux
107	integer nlev ! nombre de couches verticales
108	integer nqmax ! nombre de traceurs auxquels on applique la physique
109	integer nstep ! appel physique
110	integer julien !jour julien
111	integer itop_con(nlon)
112	integer ibas_con(nlon)
113	real gmtime
114	real pdtphys ! pas d'integration pour la physique (seconde)
115	real t_seri(nlon,nlev) ! temperature
116	real tr_seri(nlon,nlev,nbtr) ! traceur
117	real u(klon,klev)
118	real v(klon,klev)
119	real sh(nlon,nlev) ! humidite specifique
120	real rh(nlon,nlev) ! humidite relative
121	real cldliq(nlon,nlev) ! eau liquide nuageuse
122	real cldfra(nlon,nlev) ! fraction nuageuse (tous les nuages)
123	real diafra(nlon,nlev) ! fraction nuageuse (convection ou stratus artificiels)
124	real rneb(nlon,nlev) ! fraction nuageuse (grande echelle)
125	real albsol(nlon) ! albedo surface
126	real paprs(nlon,nlev+1) ! pression pour chaque inter-couche (en Pa)
127	real ps(nlon) ! pression surface
128	real pplay(nlon,nlev) ! pression pour le mileu de chaque couche (en Pa)
129	real pphi(nlon,klev) ! geopotentiel
130	real pphis(klon)
131	REAL presnivs(klev)
132	logical debutphy ! le flag de l'initialisation de la physique
133	logical lafin ! le flag de la fin de la physique
134	c Olivia
135	integer nsplit
136	REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1) !--lessivage convection
137	REAL prfl(klon,klev+1), psfl(klon,klev+1) !--lessivage large-scale
138
139	#ifdef INCA
140	REAL flxmass_w(klon,klev)
141	CHARACTER(len=8) :: solsym(nqmax)
142	#endif
143	c integer iflag_con
144
145	cAA Rem : nbtr : nombre de vrais traceurs est defini dans dimphy.h
146	c
147	c convection:
148	c -----------
149	c
150	REAL pmfu(nlon,nlev) ! flux de masse dans le panache montant
151	REAL pmfd(nlon,nlev) ! flux de masse dans le panache descendant
152	REAL pen_u(nlon,nlev) ! flux entraine dans le panache montant
153
154	c
155	c thermiques:
156	c -----------
157	c
158	real fm_therm(klon,klev+1),entr_therm(klon,klev)
159	real fm_therm1(klon,klev)
160	c
161	REAL pde_u(nlon,nlev) ! flux detraine dans le panache montant
162	REAL pen_d(nlon,nlev) ! flux entraine dans le panache descendant
163	REAL pde_d(nlon,nlev) ! flux detraine dans le panache descendant
164	c KE
165	real da(nlon,nlev),phi(nlon,nlev,nlev),mp(nlon,nlev)
166	REAL upwd(nlon,nlev) ! saturated updraft mass flux
167	REAL dnwd(nlon,nlev) ! saturated downdraft mass flux
168
169	c
170	c Couche limite:
171	c --------------
172	c
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,nqmax) ! 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	cym SAVE trs
215	REAL :: trs_tmp(klon_glo)
216
217	REAL,save,allocatable :: masktr(:,:) ! Masque reservoir de sol traceur
218	c Masque de l'echange avec la surface
219	c (1 = reservoir) ou (possible => 1 )
220	c$OMP THREADPRIVATE(masktr)
221	cym SAVE masktr
222	REAL,save,allocatable :: fshtr(:,:) ! Flux surfacique dans le reservoir de sol
223	c$OMP THREADPRIVATE(fshtr)
224	cym SAVE fshtr
225	REAL,save,allocatable :: hsoltr(:) ! Epaisseur equivalente du reservoir de sol
226	c$OMP THREADPRIVATE(hsoltr)
227	cym SAVE hsoltr
228	REAL,save,allocatable :: tautr(:) ! Constante de decroissance radioactive
229	c$OMP THREADPRIVATE(tautr)
230	cym SAVE tautr
231	REAL,save,allocatable :: vdeptr(:) ! Vitesse de depot sec dans la couche Brownienne
232	c$OMP THREADPRIVATE(vdeptr)
233	cym SAVE vdeptr
234	REAL,save,allocatable :: scavtr(:) ! Coefficient de lessivage
235	c$OMP THREADPRIVATE(scavtr)
236	cym SAVE scavtr
237	cAA
238	CHARACTER*2 itn
239	C maf ioipsl
240	CHARACTER*2 str2
241	INTEGER nhori, nvert
242	REAL zsto, zout, zjulian
243	INTEGER nid_tra
244	SAVE nid_tra
245	c$OMP THREADPRIVATE(nid_tra)
246	#ifdef INCA_AER
247	INTEGER nid_tra2,nid_tra3
248	SAVE nid_tra2,nid_tra3
249	c$OMP THREADPRIVATE(nid_tra2,nid_tra3)
250	#endif
251	c REAL x(klon,klev,nbtr+2) ! traceurs
252	INTEGER ndex(1)
253	INTEGER ndex2d(iim(jjm+1)),ndex3d(iim(jjm+1)*klev)
254	REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
255	REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
256	REAL zx_tmp_2d(iim,jjm+1), zx_tmp_3d(iim,jjm+1,klev)
257	REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1)
258	c
259	integer itau_w ! pas de temps ecriture = nstep + itau_phy
260	c
261
262	C
263	C Variables liees a l'ecriture de la bande histoire : phytrac.nc
264	c
265	c INTEGER ecrit_tra
266	c SAVE ecrit_tra
267
268	logical ok_sync
269	parameter (ok_sync = .true.)
270	C
271	C nature du traceur
272	c
273	logical,save,allocatable :: aerosol(:) ! Nature du traceur
274	c ! aerosol(it) = true => aerosol
275	c ! aerosol(it) = false => gaz
276	c ! nat_trac(it) = 1. aerosol
277	logical,save,allocatable :: clsol(:) ! clsol(it) = true => CL sol calculee
278	logical,save,allocatable :: radio(:) ! radio(it)=true => decroisssance radioactive
279	c$OMP THREADPRIVATE(aerosol,clsol,radio)
280	cym save aerosol,clsol,radio
281	C
282	c======================================================================
283	c
284	c Declaration des procedures appelees
285	c
286	c--modif convection tiedtke
287	INTEGER i, k, it
288	INTEGER iq, iiq
289	REAL delp(klon,klev)
290	c--end modif
291	c
292	c Variables liees a l'ecriture de la bande histoire physique
293	c
294	c Variables locales pour effectuer les appels en serie
295	c----------------------------------------------------
296	c
297	REAL d_tr(klon,klev), d_trs(klon) ! tendances de traceurs
298	REAL d_tr_cl(klon,klev,nbtr) ! tendance de traceurs couche limite
299	REAL d_tr_cv(klon,klev,nbtr) ! tendance de traceurs conv pour chq traceur
300	REAL d_tr_th(klon,klev,nbtr) ! la tendance des thermiques
301	REAL d_tr_dec(klon,klev,nbtr) ! la tendance de la decroissance
302	c ! radioactive du rn - > pb
303	REAL d_tr_lessi_impa(klon,klev,nbtr) ! la tendance du lessivage
304	c ! par impaction
305	REAL d_tr_lessi_nucl(klon,klev,nbtr) ! la tendance du lessivage
306	c ! par nucleation
307	REAL fluxrn(klon,klev)
308	REAL fluxpb(klon,klev)
309	REAL pbimpa(klon,klev)
310	REAL pbnucl(klon,klev)
311	REAL rn(klon,klev)
312	REAL pb(klon,klev)
313	REAL flestottr(klon,klev,nbtr) ! flux de lessivage
314	c ! dans chaque couche
315	real zmasse(klon,klev)
316	real ztra_th(klon,klev)
317
318	C
319	character*20 modname
320	character*80 abort_message
321	c
322	c Controles
323	c-------------
324	logical first,couchelimite,convection,lessivage,sorties,
325	s rnpb,inirnpb
326	save first,couchelimite,convection,lessivage,
327	s sorties,inirnpb
328	c$OMP THREADPRIVATE(first,couchelimite,convection,lessivage,
329	c$OMP+ sorties,inirnpb)
330	c data first,couchelimite,convection,lessivage,sorties
331	c s /.true.,.true.,.false.,.true.,.true./
332	c Olivia
333	data first,couchelimite,convection,lessivage,
334	s sorties
335	s /.true.,.true.,.true.,.true.,.true./
336
337
338	#ifdef INCA
339	INTEGER :: lastgas
340	INTEGER :: ncsec
341	INTEGER :: prt_flag_ts(nbtr)
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)
349	#endif
350	#ifdef INCA_AER
351	integer la
352	#endif
353	c
354	c======================================================================
355
356	#ifdef INCA
357	prt_flag_ts(:)=(/
358	#ifdef INCA_CH4
359	. 1,1,1,0,0,1,1,1,1,1,
360	. 0,1,0,0,0,0,0,1,0,0,
361	. 0,1,1,1,1,0,1,1,1,0,
362	. 1,1,1,1,1,1,1,1,1,1,
363	. 1,0,0
364	#ifdef INCA_AER
365	. ,1,1,1,1,0,1,1,1,1,0,
366	. 1,1,1,1,1,1,0,1,0,1,
367	. 1,1,1,1,0,1,0,1,1,1
368	#endif
369	#endif
370	#ifdef INCA_NMHC
371	. 1,1,1,1,1,1,1,1,1,1,
372	. 1,1,1,1,1,1,1,1,1,1,
373	. 1,1,1,1,1,1,1,1,1,1,
374	. 1,1,1,1,1,1,1,1,1,1,
375	. 1,1,1,1,1,1,1,1,1,1,
376	. 1,1,1,1,1,1,1,1,1,1,
377	. 1,1,1,1,1,1,1,1,1,1,
378	. 1,1,1,1,1,1,1,1,1,1,
379	. 1,1,1,1,1,1,1
380	#ifdef INCA_AER
381	. ,1,1,1,1,0,1,1,1,1,0,
382	. 1,1,1,1,1,1,0,1,0,1,
383	. 1,1,1,1,0,1,0,1,1,1
384	#endif
385	#endif
386	#if defined(INCA_AER) && !defined(INCA_CH4) && !defined(INCA_NMHC)
387	. 1,1,1,1,1,1,1,1,1,1,
388	. 1,1,1,1,1,1,1,1,1,1,
389	. 1,1,1,1,1,1,1,1,1
390	#endif
391	#if defined(INCA) && !defined(INCA_CH4) && !defined(INCA_NMHC) && !defined(INCA_AER)
392	. 1,1,1,1,1,1,1,1,1,1,
393	. 1
394	#endif
395
396	. /)
397	dummy(:,:) = 0.
398
399	#endif
400	modname='phytrac'
401
402	ps(:)=paprs(:,1)
403
404	if (debutphy) then
405	allocate( trs(klon,nbtr) )
406	allocate( masktr(klon,nbtr))
407	allocate( fshtr(klon,nbtr) )
408	allocate( hsoltr(nbtr))
409	allocate( tautr(nbtr))
410	allocate( vdeptr(nbtr))
411	allocate( scavtr(nbtr))
412	allocate( aerosol(nbtr))
413	allocate( clsol(nbtr))
414	allocate( radio(nbtr))
415
416
417	ecrit_tra = NINT(86400./pdtphys *ecritphy)
418	print*,'dans phytrac ',pdtphys,ecritphy,ecrit_tra
419
420	if(nbtr.lt.nqmax) then
421	c print*,'NQMAX=',nqmax
422	c print*,'NBTR=',nbtr
423	abort_message='See above'
424	call abort_gcm(modname,abort_message,1)
425	endif
426
427	inirnpb=rnpb
428	PRINT*, 'La frequence de sortie traceurs est ', ecrit_tra
429	C
430	c=============================================================
431	c Initialisation des sorties
432	c=============================================================
433
434	#ifdef CPP_IOIPSL
435	#include "ini_histrac.h"
436	#endif
437
438	c======================================================================
439	c Initialisation de certaines variables pour le Rn et le Pb
440	c======================================================================
441
442	c Initialisation du traceur dans le sol (couche limite radonique)
443	c
444	c print*,'valeur de debut dans phytrac :',debutphy
445	trs(:,:) = 0.
446	c$OMP MASTER
447	if (is_mpi_root) then
448	trs_tmp(:)=0.
449	open (99,file='starttrac',status='old',
450	. err=999,form='formatted')
451	read(99,*) (trs_tmp(i),i=1,klon_glo)
452	999 close(99)
453	endif
454	c$OMP END MASTER
455	call Scatter(trs_tmp,trs(:,1))
456
457	c print*, 'apres starttrac'
458
459	c Initialisation de la fraction d'aerosols lessivee
460	c
461	d_tr_lessi_impa(:,:,:) = 0.
462	d_tr_lessi_nucl(:,:,:) = 0.
463	c
464	c Initialisation de la nature des traceurs
465	c
466	DO it = 1, nqmax
467	aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut
468	radio(it) = .FALSE. ! Par defaut pas de passage par radiornpb
469	clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit
470	ENDDO
471	c
472	ENDIF ! fin debutphy
473	c Initialisation du traceur dans le sol (couche limite radonique)
474	if(inirnpb) THEN
475	c
476	radio(1)= .true.
477	radio(2)= .true.
478	clsol(1)= .true.
479	clsol(2)= .true.
480	aerosol(2) = .TRUE. ! le Pb est un aerosol
481	c
482	call initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr
483	. ,vdeptr,scavtr)
484	inirnpb=.false.
485	endif
486	#ifdef INCA
487	call VTe(VTphysiq)
488	call VTb(VTinca)
489	!======================================================================
490	! Chimie
491	!======================================================================
492
493	calday = FLOAT(julien) + gmtime
494	ncsec = NINT (86400.*gmtime)
495
496	DO k = 1, nlev
497	pdel(:,k) = paprs(:,k) - paprs (:,k+1)
498	END DO
499
500	#ifdef INCAINFO
501	PRINT *, 'CHEMMAIN @ ', calday, ' ... '
502	DO it = 1, nbtr
503	PRINT *, solsym(it), MINVAL(tr_seri(:,:,it)),
504	$ MAXVAL(tr_seri(:,:,it))
505	END DO
506	#endif
507
508
509	#ifdef INCA_AER
510
511	! Changement Anne 01/04/2005
512	CALL aerosolmain (tr_seri,
513	$ pdtphys,
514	$ pplay,
515	$ pdel,
516	$ prfl,
517	$ pmflxr,
518	$ psfl,
519	$ pmflxs,
520	$ pmfu,
521	$ itop_con,
522	$ ibas_con,
523	$ pphi,
524	$ airephy, ! paire,
525	$ nstep,
526	$ rneb, ! for chimiaq
527	$ t_seri, ! for chimiaq
528	$ rh,
529	$ lafin)
530	! fin changement anne
531
532	#endif
533
534	CALL chemmain (tr_seri, !mmr
535	$ nstep, !nstep
536	$ calday, !calday
537	$ julien, !ncdate
538	$ ncsec, !ncsec
539	$ 1, !lat
540	$ pdtphys, !delt
541	$ paprs(1,1), !ps
542	$ pplay, !pmid
543	$ pdel, !pdel
544	$ airephy,
545	$ pctsrf(1,1),!oro
546	$ ftsol, !tsurf
547	$ albsol, !albs
548	$ pphi, !zma
549	$ pphis, !phis
550	$ cldfra, !cldfr
551	$ rneb, !cldfr_st
552	$ diafra, !cldfr_cv
553	$ itop_con, !cldtop
554	$ ibas_con, !cldbot
555	$ cldliq, !cwat
556	$ prfl, !flxrst
557	$ pmflxr, !flxrcv
558	$ psfl, !flxsst
559	$ pmflxs, !flxscv
560	$ pmfu, !flxupd
561	$ flxmass_w, !flxmass_w
562	$ t_seri, !tfld
563	$ sh, !sh
564	$ rh, !rh
565	$ .false., !wrhstts
566	$ hbuf, !hbuf
567	$ obuf, !obuf
568	$ iip1, !nx
569	$ jjp1, !ny
570	$ source,
571	$ solsym)
572	#ifdef INCAINFO
573	#ifdef INCA_AER
574
575	c Budget calculation for aerosol species
576	CALL tbudget(airephy,pdtphys,nstep,tr_seri,.false.)
577
578	c-- summary info----------------------------------------------------------------
579
580	if (MOD(nstep,nint(86400./pdtphys)) .eq. 0) then
581	print *, "global aerosol optical thickness "
582
583	write (form,'(A,I2,A)') "(A,",trnx-trmx+1,"A10)"
584	print form,"lamba [nm] ", (solsym(it),it=trmx,trnx)
585
586	write (form,'(A,I2,A)') "(I11,",trnx-trmx+1,"F10.4)"
587	do i=1,las
588	print form,int(lambda(i)),(sum(tausum(:,i,it)*airephy)/sum(airephy),it=trmx,trnx)
589	enddo
590
591	print ,"global mean angstroem component ", sum(angstairephy)/sum(airephy)
592	endif
593	#endif
594	#endif
595
596	#ifdef INCAINFO
597	PRINT *, 'OK.'
598	DO it = 1, nbtr
599	PRINT *, solsym(it), MINVAL(tr_seri(:,:,it)),
600	$ MAXVAL(tr_seri(:,:,it))
601	END DO
602	#endif
603	call VTe(VTinca)
604	call VTb(VTphysiq)
605	#else
606
607	c Abder
608	ctestmaf if(nqmax.gt.2) aerosol(3)=.true.
609
610	do i=1,nlon
611	pftsol1(i) = ftsol(i,1)
612	pftsol2(i) = ftsol(i,2)
613	pftsol3(i) = ftsol(i,3)
614	pftsol4(i) = ftsol(i,4)
615
616	ppsrf1(i) = pctsrf(i,1)
617	ppsrf2(i) = pctsrf(i,2)
618	ppsrf3(i) = pctsrf(i,3)
619	ppsrf4(i) = pctsrf(i,4)
620
621	enddo
622	c Abder
623	#endif
624	c======================================================================
625	c Calcul de l'effet de la convection
626	c======================================================================
627	c print*,'Avant convection'
628	c do it=1,nqmax
629	c WRITE(itn,'(i2)') it
630	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
631	c enddo
632
633	if (convection) then
634
635	c print*,'Pas de temps dans phytrac : ',pdtphys
636	DO it=1, nqmax
637	#ifdef INCA
638	IF ( conv_flg(it) == 0 ) CYCLE
639	#endif
640	if (iflag_con.eq.2) then
641	c tiedke
642	CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
643	. pplay, paprs, tr_seri(1,1,it), d_tr_cv(1,1,it))
644	else if (iflag_con.eq.3) then
645	c KE
646	call cvltr(pdtphys, da, phi, mp, paprs,pplay, tr_seri(1,1,it),
647	. upwd,dnwd,d_tr_cv(1,1,it))
648	endif
649
650	DO k = 1, nlev
651	DO i = 1, klon
652	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it)
653	ENDDO
654	ENDDO
655	#ifdef INCA
656	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '
657	. //solsym(it))
658	#else
659	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '//itn)
660	#endif
661	ENDDO
662	c print*,'apres nflxtr'
663
664	endif ! convection
665	c print*,'Apres convection'
666	c do it=1,nqmax
667	c WRITE(itn,'(i1)') it
668	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
669	c enddo
670
671
672	c======================================================================
673	c Calcul de l'effet des thermiques
674	c======================================================================
675
676	do k=1,klev
677	do i=1,klon
678	zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
679	enddo
680	enddo
681
682	c print*,'masse dans ph ',zmasse
683	do it=1,nqmax
684	do k=1,klev
685	do i=1,klon
686	d_tr_th(i,k,it)=0.
687	tr_seri(i,k,it)=max(tr_seri(i,k,it),0.)
688	tr_seri(i,k,it)=min(tr_seri(i,k,it),1.e10)
689	enddo
690	enddo
691	enddo
692
693	if (iflag_thermals.gt.0) then
694	c print*,'calcul de leffet des thermiques'
695	nsplit=10
696	DO it=1, nqmax
697	c WRITE(itn,'(i1)') it
698	c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'conv it='//itn)
699	c print*,'avant dqthermiquesretro'
700	c call dump2d(iim,jjm-1,tr_seri(2,1,1),'TR_SERI ')
701
702	do isplit=1,nsplit
703	c Abderr 25 11 02
704	C Thermiques
705	c print*,'Avant dans phytrac'
706	call dqthermcell(klon,klev,pdtphys/nsplit
707	. ,fm_therm,entr_therm,zmasse
708	. ,tr_seri(1:klon,1:klev,it),d_tr,ztra_th)
709
710	do k=1,klev
711	do i=1,klon
712	d_tr(i,k)=pdtphys*d_tr(i,k)/nsplit
713	d_tr_th(i,k,it)=d_tr_th(i,k,it)+d_tr(i,k)
714	tr_seri(i,k,it)=max(tr_seri(i,k,it)+d_tr(i,k),0.)
715	enddo
716	enddo
717	enddo ! nsplit
718	c print*,'apres thermiques'
719	c call dump2d(iim,jjm-1,d_tr_th(1,1,1),'d_tr_th ')
720	c do k=1,klev
721	c print*,'d_tr_th(',k,')=',tr_seri(280,k,1)
722	c enddo
723
724	c WRITE(itn,'(i1)') it
725	c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'therm it='//itn)
726	ENDDO ! it
727	endif ! Thermiques
728	c print*,'ATTENTION: sdans thermniques'
729
730	c======================================================================
731	c Calcul de l'effet de la couche limite
732	c======================================================================
733	c print *,'Avant couchelimite'
734	c do it=1,nqmax
735	c WRITE(itn,'(i1)') it
736	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
737	c enddo
738
739	if (couchelimite) then
740
741	DO k = 1, nlev
742	DO i = 1, klon
743	delp(i,k) = paprs(i,k)-paprs(i,k+1)
744	ENDDO
745	ENDDO
746
747	C maf modif pour tenir compte du cas rnpb + traceur
748	DO it=1, nqmax
749	#ifdef INCA
750	IF ( pbl_flg(it) == 0 ) CYCLE
751	#endif
752	c print *,'it',it,clsol(it)
753	if (clsol(it)) then ! couche limite avec quantite dans le sol calculee
754	CALL cltracrn(it, pdtphys, yu1, yv1,
755	e coefh,t_seri,ftsol,pctsrf,
756	e tr_seri(1,1,it),trs(1,it),
757	e paprs, pplay, delp,
758	e masktr(1,it),fshtr(1,it),hsoltr(it),
759	e tautr(it),vdeptr(it),
760	e xlat,
761	s d_tr_cl(1,1,it),d_trs)
762	DO k = 1, nlev
763	DO i = 1, klon
764	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
765	ENDDO
766	ENDDO
767	c
768	c Traceur ds sol
769	c
770	DO i = 1, klon
771	trs(i,it) = trs(i,it) + d_trs(i)
772	END DO
773	C
774	C maf provisoire suppression des prints
775	C WRITE(itn,'(i1)') it
776	C CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
777	else ! couche limite avec flux prescrit
778	#ifndef INCA
779
780	Cmaf provisoire source / traceur a creer
781	DO i=1, klon
782	source(i,it) = 0.0 ! pas de source, pour l'instant
783	ENDDO
784	C
785	#endif
786	CALL cltrac(pdtphys, coefh,t_seri,
787	s tr_seri(1,1,it), source(:,it),
788	e paprs, pplay, delp,
789	s d_tr_cl(1,1,it))
790	DO k = 1, nlev
791	DO i = 1, klon
792	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
793	ENDDO
794	ENDDO
795	Cmaf WRITE(itn,'(i1)') it
796	cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
797	endif
798	ENDDO
799	c
800	endif ! couche limite
801
802	c print*,'Apres couchelimite'
803	c do it=1,nqmax
804	c WRITE(itn,'(i1)') it
805	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
806	c enddo
807
808	c======================================================================
809	c Calcul de l'effet du puits radioactif
810	c======================================================================
811
812	C MAF il faudrait faire une modification pour passer dans radiornpb
813	c si radio=true mais pour l'instant radiornpb propre au cas rnpb
814	if(rnpb) then
815	print *, 'decroissance radiactive activee'
816	call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
817	C
818	DO it=1,nqmax
819	if(radio(it)) then
820	DO k = 1, nlev
821	DO i = 1, klon
822	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
823	ENDDO
824	ENDDO
825	WRITE(itn,'(i1)') it
826	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn)
827	endif
828	ENDDO
829	c
830	endif ! rnpb decroissance radioactive
831	C
832	c======================================================================
833	c Calcul de l'effet de la precipitation
834	c======================================================================
835
836	c print*,'LESSIVAGE =',lessivage
837	IF (lessivage) THEN
838
839	c print*,'avant lessivage'
840
841	d_tr_lessi_nucl(:,:,:) = 0.
842	d_tr_lessi_impa(:,:,:) = 0.
843	flestottr(:,:,:) = 0.
844	c
845	c tendance des aerosols nuclees et impactes
846	c
847	DO it = 1, nqmax
848	IF (aerosol(it)) THEN
849	DO k = 1, nlev
850	DO i = 1, klon
851	d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +
852	s ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
853	d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +
854	s ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
855	ENDDO
856	ENDDO
857	ENDIF
858	ENDDO
859	c
860	c Mises a jour des traceurs + calcul des flux de lessivage
861	c Mise a jour due a l'impaction et a la nucleation
862	c
863	c call dump2d(iim,jjm-1,frac_impa(2:klon-1,10),'FRACIMPA')
864	c call dump2d(iim,jjm-1,frac_nucl(2:klon-1,10),'FRACNUCL')
865	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3')
866	DO it = 1, nqmax
867	c print*,'IT=',it,aerosol(it)
868	IF (aerosol(it)) THEN
869	c print*,'IT=',it,' On lessive'
870	DO k = 1, nlev
871	DO i = 1, klon
872	tr_seri(i,k,it)=tr_seri(i,k,it)
873	s frac_impa(i,k)frac_nucl(i,k)
874	ENDDO
875	ENDDO
876	ENDIF
877	ENDDO
878	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3B')
879	c
880	c Flux lessivage total
881	c
882	DO it = 1, nqmax
883	DO k = 1, nlev
884	DO i = 1, klon
885	flestottr(i,k,it) = flestottr(i,k,it) -
886	s ( d_tr_lessi_nucl(i,k,it) +
887	s d_tr_lessi_impa(i,k,it) ) *
888	s ( paprs(i,k)-paprs(i,k+1) ) /
889	s (RG * pdtphys)
890	ENDDO
891	ENDDO
892	c
893	Cmaf WRITE(itn,'(i1)') it
894	Cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
895	ENDDO
896	c
897	c print*,'apres lessivage'
898	ENDIF
899	Cc
900	DO k = 1, nlev
901	DO i = 1, klon
902	fluxrn(i,k) = flestottr(i,k,1)
903	fluxpb(i,k) = flestottr(i,k,2)
904	rn(i,k) = tr_seri(i,k,1)
905	pb(i,k) = tr_seri(i,k,2)
906	pbnucl(i,k)=d_tr_lessi_nucl(i,k,2)
907	pbimpa(i,k)=d_tr_lessi_impa(i,k,2)
908	ENDDO
909	ENDDO
910
911	c=============================================================
912	c Ecriture des sorties
913	c=============================================================
914
915	#ifdef CPP_IOIPSL
916	#include "write_histrac.h"
917	#endif
918
919	c=============================================================
920
921	if (lafin) then
922	print*, 'c est la fin de la physique'
923	call Gather(trs(:,1),trs_tmp)
924	c$OMP MASTER
925	if (is_mpi_root) then
926
927	open (99,file='restarttrac', form='formatted')
928	do i=1,klon_glo
929	write(99,*) trs_tmp(i)
930	enddo
931	PRINT*, 'Ecriture du fichier restarttrac'
932	close(99)
933	endif
934	c$OMP END MASTER
935	else
936	c print*, 'physique pas fini'
937	endif
938
939
940	RETURN
941	END

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

Download in other formats:

Original Format