Context Navigation

phytrac.F @ 901

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

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

Download in other formats:

Original Format