Context Navigation

phytrac.F @ 1153

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

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

Download in other formats:

Original Format