Context Navigation

phytrac.F @ 715

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

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

Download in other formats:

Original Format