Context Navigation

phytrac.F @ 1161

Last change on this file since 1161 was 1150, checked in by jghattas, 15 years ago

Ajoute possiblilite de forcer avec d'autre aersols avec le meme principe que pour les so4.

Par default la convergence est rempu avec flag_aersosol=1 + ok_ade=ok_aie=y. Exactement les memes resultats peuvent etre retrouver avec new_aod=n. Les resultats sont exactement les memes sans prise en compte aucun aerosol, avec ok_ade=ok_aie=n.

flag_aerosol indique quel aersosol ou combinasion a utiliser (=1 uniquement les SO4 comme avant)
les fichiers d'entree so4.runXXXX.cdf change du nom pour majescule SO4.runXXXX.cdf

readsulfate change du nom pour readaerosol qui trait tous les aerosols
radlwsw change du nom pour radlwsw_aero.
aeropt_2bands.F90 et aeropt_5wv.F90 correspond a un reecriture de aeropt.F (premier et deuxieme moitie respectivement)
aeropt.F est gardé pour retrouver la convergence si demande avec new_aod=false (=true default)
sw_aero est un version evolue de sw_lmdar4 (dans fichier radiation_ar4.F)

ACo + JG

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 26.4 KB

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

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

Download in other formats:

Original Format