Context Navigation

phytrac.F @ 815

Last change on this file since 815 was 814, checked in by Laurent Fairhead, 17 years ago
Rajout de la physique utilisant les thermiques FH 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	c jg: c'est ca qu'on veut?????
418	ecrit_tra = FLOAT(NINT(86400./pdtphys *ecritphy))
419	print*,'dans phytrac ',pdtphys,ecritphy,ecrit_tra
420
421	if(nbtr.lt.nqmax) then
422	c print*,'NQMAX=',nqmax
423	c print*,'NBTR=',nbtr
424	abort_message='See above'
425	call abort_gcm(modname,abort_message,1)
426	endif
427
428	inirnpb=rnpb
429	PRINT*, 'La frequence de sortie traceurs est ', ecrit_tra
430	C
431	c=============================================================
432	c Initialisation des sorties
433	c=============================================================
434
435	#ifdef CPP_IOIPSL
436	#include "ini_histrac.h"
437	#endif
438
439	c======================================================================
440	c Initialisation de certaines variables pour le Rn et le Pb
441	c======================================================================
442
443	c Initialisation du traceur dans le sol (couche limite radonique)
444	c
445	c print*,'valeur de debut dans phytrac :',debutphy
446	trs(:,:) = 0.
447	c$OMP MASTER
448	if (is_mpi_root) then
449	trs_tmp(:)=0.
450	open (99,file='starttrac',status='old',
451	. err=999,form='formatted')
452	read(99,*) (trs_tmp(i),i=1,klon_glo)
453	999 close(99)
454	endif
455	c$OMP END MASTER
456	call Scatter(trs_tmp,trs(:,1))
457
458	c print*, 'apres starttrac'
459
460	c Initialisation de la fraction d'aerosols lessivee
461	c
462	d_tr_lessi_impa(:,:,:) = 0.
463	d_tr_lessi_nucl(:,:,:) = 0.
464	c
465	c Initialisation de la nature des traceurs
466	c
467	DO it = 1, nqmax
468	aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut
469	radio(it) = .FALSE. ! Par defaut pas de passage par radiornpb
470	clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit
471	ENDDO
472	c
473	ENDIF ! fin debutphy
474	c Initialisation du traceur dans le sol (couche limite radonique)
475	if(inirnpb) THEN
476	c
477	radio(1)= .true.
478	radio(2)= .true.
479	clsol(1)= .true.
480	clsol(2)= .true.
481	aerosol(2) = .TRUE. ! le Pb est un aerosol
482	c
483	call initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr
484	. ,vdeptr,scavtr)
485	inirnpb=.false.
486	endif
487	#ifdef INCA
488	call VTe(VTphysiq)
489	call VTb(VTinca)
490	!======================================================================
491	! Chimie
492	!======================================================================
493
494	calday = FLOAT(julien) + gmtime
495	ncsec = NINT (86400.*gmtime)
496
497	DO k = 1, nlev
498	pdel(:,k) = paprs(:,k) - paprs (:,k+1)
499	END DO
500
501	#ifdef INCAINFO
502	PRINT *, 'CHEMMAIN @ ', calday, ' ... '
503	DO it = 1, nbtr
504	PRINT *, solsym(it), MINVAL(tr_seri(:,:,it)),
505	$ MAXVAL(tr_seri(:,:,it))
506	END DO
507	#endif
508
509
510	#ifdef INCA_AER
511
512	! Changement Anne 01/04/2005
513	CALL aerosolmain (tr_seri,
514	$ pdtphys,
515	$ pplay,
516	$ pdel,
517	$ prfl,
518	$ pmflxr,
519	$ psfl,
520	$ pmflxs,
521	$ pmfu,
522	$ itop_con,
523	$ ibas_con,
524	$ pphi,
525	$ airephy, ! paire,
526	$ nstep,
527	$ rneb, ! for chimiaq
528	$ t_seri, ! for chimiaq
529	$ rh,
530	$ lafin)
531	! fin changement anne
532
533	#endif
534
535	CALL chemmain (tr_seri, !mmr
536	$ nstep, !nstep
537	$ calday, !calday
538	$ julien, !ncdate
539	$ ncsec, !ncsec
540	$ 1, !lat
541	$ pdtphys, !delt
542	$ paprs(1,1), !ps
543	$ pplay, !pmid
544	$ pdel, !pdel
545	$ airephy,
546	$ pctsrf(1,1),!oro
547	$ ftsol, !tsurf
548	$ albsol, !albs
549	$ pphi, !zma
550	$ pphis, !phis
551	$ cldfra, !cldfr
552	$ rneb, !cldfr_st
553	$ diafra, !cldfr_cv
554	$ itop_con, !cldtop
555	$ ibas_con, !cldbot
556	$ cldliq, !cwat
557	$ prfl, !flxrst
558	$ pmflxr, !flxrcv
559	$ psfl, !flxsst
560	$ pmflxs, !flxscv
561	$ pmfu, !flxupd
562	$ flxmass_w, !flxmass_w
563	$ t_seri, !tfld
564	$ sh, !sh
565	$ rh, !rh
566	$ .false., !wrhstts
567	$ hbuf, !hbuf
568	$ obuf, !obuf
569	$ iip1, !nx
570	$ jjp1, !ny
571	$ source,
572	$ solsym)
573	#ifdef INCAINFO
574	#ifdef INCA_AER
575
576	c Budget calculation for aerosol species
577	CALL tbudget(airephy,pdtphys,nstep,tr_seri,.false.)
578
579	c-- summary info----------------------------------------------------------------
580
581	if (MOD(nstep,nint(86400./pdtphys)) .eq. 0) then
582	print *, "global aerosol optical thickness "
583
584	write (form,'(A,I2,A)') "(A,",trnx-trmx+1,"A10)"
585	print form,"lamba [nm] ", (solsym(it),it=trmx,trnx)
586
587	write (form,'(A,I2,A)') "(I11,",trnx-trmx+1,"F10.4)"
588	do i=1,las
589	print form,int(lambda(i)),(sum(tausum(:,i,it)*airephy)/sum(airephy),it=trmx,trnx)
590	enddo
591
592	print ,"global mean angstroem component ", sum(angstairephy)/sum(airephy)
593	endif
594	#endif
595	#endif
596
597	#ifdef INCAINFO
598	PRINT *, 'OK.'
599	DO it = 1, nbtr
600	PRINT *, solsym(it), MINVAL(tr_seri(:,:,it)),
601	$ MAXVAL(tr_seri(:,:,it))
602	END DO
603	#endif
604	call VTe(VTinca)
605	call VTb(VTphysiq)
606	#else
607
608	c Abder
609	ctestmaf if(nqmax.gt.2) aerosol(3)=.true.
610
611	do i=1,nlon
612	pftsol1(i) = ftsol(i,1)
613	pftsol2(i) = ftsol(i,2)
614	pftsol3(i) = ftsol(i,3)
615	pftsol4(i) = ftsol(i,4)
616
617	ppsrf1(i) = pctsrf(i,1)
618	ppsrf2(i) = pctsrf(i,2)
619	ppsrf3(i) = pctsrf(i,3)
620	ppsrf4(i) = pctsrf(i,4)
621
622	enddo
623	c Abder
624	#endif
625	c======================================================================
626	c Calcul de l'effet de la convection
627	c======================================================================
628	c print*,'Avant convection'
629	c do it=1,nqmax
630	c WRITE(itn,'(i2)') it
631	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
632	c enddo
633
634	if (convection) then
635
636	c print*,'Pas de temps dans phytrac : ',pdtphys
637	DO it=1, nqmax
638	#ifdef INCA
639	IF ( conv_flg(it) == 0 ) CYCLE
640	#endif
641	if (iflag_con.eq.2) then
642	c tiedke
643	CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
644	. pplay, paprs, tr_seri(1,1,it), d_tr_cv(1,1,it))
645	else if (iflag_con.eq.3) then
646	c KE
647	call cvltr(pdtphys, da, phi, mp, paprs,pplay, tr_seri(1,1,it),
648	. upwd,dnwd,d_tr_cv(1,1,it))
649	endif
650
651	DO k = 1, nlev
652	DO i = 1, klon
653	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it)
654	ENDDO
655	ENDDO
656	#ifdef INCA
657	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '
658	. //solsym(it))
659	#else
660	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '//itn)
661	#endif
662	ENDDO
663	c print*,'apres nflxtr'
664
665	endif ! convection
666	c print*,'Apres convection'
667	c do it=1,nqmax
668	c WRITE(itn,'(i1)') it
669	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
670	c enddo
671
672
673	c======================================================================
674	c Calcul de l'effet des thermiques
675	c======================================================================
676
677	do k=1,klev
678	do i=1,klon
679	zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
680	enddo
681	enddo
682
683	c print*,'masse dans ph ',zmasse
684	do it=1,nqmax
685	do k=1,klev
686	do i=1,klon
687	d_tr_th(i,k,it)=0.
688	tr_seri(i,k,it)=max(tr_seri(i,k,it),0.)
689	tr_seri(i,k,it)=min(tr_seri(i,k,it),1.e10)
690	enddo
691	enddo
692	enddo
693
694	if (iflag_thermals.gt.0) then
695	c print*,'calcul de leffet des thermiques'
696	nsplit=10
697	DO it=1, nqmax
698	c WRITE(itn,'(i1)') it
699	c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'conv it='//itn)
700	c print*,'avant dqthermiquesretro'
701	c call dump2d(iim,jjm-1,tr_seri(2,1,1),'TR_SERI ')
702
703	do isplit=1,nsplit
704	c Abderr 25 11 02
705	C Thermiques
706	c print*,'Avant dans phytrac'
707	call dqthermcell(klon,klev,pdtphys/nsplit
708	. ,fm_therm,entr_therm,zmasse
709	. ,tr_seri(1:klon,1:klev,it),d_tr,ztra_th)
710
711	do k=1,klev
712	do i=1,klon
713	d_tr(i,k)=pdtphys*d_tr(i,k)/nsplit
714	d_tr_th(i,k,it)=d_tr_th(i,k,it)+d_tr(i,k)
715	tr_seri(i,k,it)=max(tr_seri(i,k,it)+d_tr(i,k),0.)
716	enddo
717	enddo
718	enddo ! nsplit
719	c print*,'apres thermiques'
720	c call dump2d(iim,jjm-1,d_tr_th(1,1,1),'d_tr_th ')
721	c do k=1,klev
722	c print*,'d_tr_th(',k,')=',tr_seri(280,k,1)
723	c enddo
724
725	c WRITE(itn,'(i1)') it
726	c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'therm it='//itn)
727	ENDDO ! it
728	endif ! Thermiques
729	c print*,'ATTENTION: sdans thermniques'
730
731	c======================================================================
732	c Calcul de l'effet de la couche limite
733	c======================================================================
734	c print *,'Avant couchelimite'
735	c do it=1,nqmax
736	c WRITE(itn,'(i1)') it
737	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
738	c enddo
739
740	if (couchelimite) then
741
742	DO k = 1, nlev
743	DO i = 1, klon
744	delp(i,k) = paprs(i,k)-paprs(i,k+1)
745	ENDDO
746	ENDDO
747
748	C maf modif pour tenir compte du cas rnpb + traceur
749	DO it=1, nqmax
750	#ifdef INCA
751	IF ( pbl_flg(it) == 0 ) CYCLE
752	#endif
753	c print *,'it',it,clsol(it)
754	if (clsol(it)) then ! couche limite avec quantite dans le sol calculee
755	CALL cltracrn(it, pdtphys, yu1, yv1,
756	e coefh,t_seri,ftsol,pctsrf,
757	e tr_seri(1,1,it),trs(1,it),
758	e paprs, pplay, delp,
759	e masktr(1,it),fshtr(1,it),hsoltr(it),
760	e tautr(it),vdeptr(it),
761	e xlat,
762	s d_tr_cl(1,1,it),d_trs)
763	DO k = 1, nlev
764	DO i = 1, klon
765	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
766	ENDDO
767	ENDDO
768	c
769	c Traceur ds sol
770	c
771	DO i = 1, klon
772	trs(i,it) = trs(i,it) + d_trs(i)
773	END DO
774	C
775	C maf provisoire suppression des prints
776	C WRITE(itn,'(i1)') it
777	C CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
778	else ! couche limite avec flux prescrit
779	#ifndef INCA
780
781	Cmaf provisoire source / traceur a creer
782	DO i=1, klon
783	source(i,it) = 0.0 ! pas de source, pour l'instant
784	ENDDO
785	C
786	#endif
787	CALL cltrac(pdtphys, coefh,t_seri,
788	s tr_seri(1,1,it), source(:,it),
789	e paprs, pplay, delp,
790	s d_tr_cl(1,1,it))
791	DO k = 1, nlev
792	DO i = 1, klon
793	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
794	ENDDO
795	ENDDO
796	Cmaf WRITE(itn,'(i1)') it
797	cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
798	endif
799	ENDDO
800	c
801	endif ! couche limite
802
803	c print*,'Apres couchelimite'
804	c do it=1,nqmax
805	c WRITE(itn,'(i1)') it
806	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
807	c enddo
808
809	c======================================================================
810	c Calcul de l'effet du puits radioactif
811	c======================================================================
812
813	C MAF il faudrait faire une modification pour passer dans radiornpb
814	c si radio=true mais pour l'instant radiornpb propre au cas rnpb
815	if(rnpb) then
816	c print *, 'decroissance radiactive activee'
817	call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
818	C
819	DO it=1,nqmax
820	if(radio(it)) then
821	DO k = 1, nlev
822	DO i = 1, klon
823	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
824	ENDDO
825	ENDDO
826	WRITE(itn,'(i1)') it
827	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn)
828	endif
829	ENDDO
830	c
831	endif ! rnpb decroissance radioactive
832	C
833	c======================================================================
834	c Calcul de l'effet de la precipitation
835	c======================================================================
836
837	c print*,'LESSIVAGE =',lessivage
838	IF (lessivage) THEN
839
840	c print*,'avant lessivage'
841
842	d_tr_lessi_nucl(:,:,:) = 0.
843	d_tr_lessi_impa(:,:,:) = 0.
844	flestottr(:,:,:) = 0.
845	c
846	c tendance des aerosols nuclees et impactes
847	c
848	DO it = 1, nqmax
849	IF (aerosol(it)) THEN
850	DO k = 1, nlev
851	DO i = 1, klon
852	d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +
853	s ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
854	d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +
855	s ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
856	ENDDO
857	ENDDO
858	ENDIF
859	ENDDO
860	c
861	c Mises a jour des traceurs + calcul des flux de lessivage
862	c Mise a jour due a l'impaction et a la nucleation
863	c
864	c call dump2d(iim,jjm-1,frac_impa(2:klon-1,10),'FRACIMPA')
865	c call dump2d(iim,jjm-1,frac_nucl(2:klon-1,10),'FRACNUCL')
866	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3')
867	DO it = 1, nqmax
868	c print*,'IT=',it,aerosol(it)
869	IF (aerosol(it)) THEN
870	c print*,'IT=',it,' On lessive'
871	DO k = 1, nlev
872	DO i = 1, klon
873	tr_seri(i,k,it)=tr_seri(i,k,it)
874	s frac_impa(i,k)frac_nucl(i,k)
875	ENDDO
876	ENDDO
877	ENDIF
878	ENDDO
879	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3B')
880	c
881	c Flux lessivage total
882	c
883	DO it = 1, nqmax
884	DO k = 1, nlev
885	DO i = 1, klon
886	flestottr(i,k,it) = flestottr(i,k,it) -
887	s ( d_tr_lessi_nucl(i,k,it) +
888	s d_tr_lessi_impa(i,k,it) ) *
889	s ( paprs(i,k)-paprs(i,k+1) ) /
890	s (RG * pdtphys)
891	ENDDO
892	ENDDO
893	c
894	Cmaf WRITE(itn,'(i1)') it
895	Cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
896	ENDDO
897	c
898	c print*,'apres lessivage'
899	ENDIF
900	Cc
901	DO k = 1, nlev
902	DO i = 1, klon
903	fluxrn(i,k) = flestottr(i,k,1)
904	fluxpb(i,k) = flestottr(i,k,2)
905	rn(i,k) = tr_seri(i,k,1)
906	pb(i,k) = tr_seri(i,k,2)
907	pbnucl(i,k)=d_tr_lessi_nucl(i,k,2)
908	pbimpa(i,k)=d_tr_lessi_impa(i,k,2)
909	ENDDO
910	ENDDO
911
912	c=============================================================
913	c Ecriture des sorties
914	c=============================================================
915
916	#ifdef CPP_IOIPSL
917	#include "write_histrac.h"
918	#endif
919
920	c=============================================================
921
922	if (lafin) then
923	print*, 'c est la fin de la physique'
924	call Gather(trs(:,1),trs_tmp)
925	c$OMP MASTER
926	if (is_mpi_root) then
927
928	open (99,file='restarttrac', form='formatted')
929	do i=1,klon_glo
930	write(99,*) trs_tmp(i)
931	enddo
932	PRINT*, 'Ecriture du fichier restarttrac'
933	close(99)
934	endif
935	c$OMP END MASTER
936	else
937	c print*, 'physique pas fini'
938	endif
939
940
941	RETURN
942	END

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

Download in other formats:

Original Format