Context Navigation

phytrac.F @ 634

Last change on this file since 634 was 634, checked in by Laurent Fairhead, 19 years ago
Modifications faites à la physique pour la rendre parallele YM Une branche de travail LMDZ4_par_0 a été créée provisoirement afin de tester les modifs pleinement avant leurs inclusions dans le tronc principal LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 28.7 KB

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

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

Download in other formats:

Original Format