Context Navigation

phytrac.F @ 95

Last change on this file since 95 was 79, checked in by (none), 25 years ago
This commit was manufactured by cvs2svn to create branch 'rel-LF'.
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 17.5 KB

Rev	Line
[2]	1	SUBROUTINE phytrac (rnpb,
	2	I debutphy,
	3	I nqmax,
	4	I nlon,nlev,pdtphys,
	5	I t_seri,paprs,pplay,
	6	I pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
	7	I coefh,yu1,yv1,ftsol,pctsrf,xlat,
	8	I frac_impa, frac_nucl,
	9	I xlon,presnivs,paire,pphis,
	10	O tr_seri)
	11	USE ioipsl
	12	IMPLICIT none
	13	c======================================================================
	14	c Auteur(s) FH
	15	c Objet: Moniteur general des tendances traceurs
	16	c
	17	cAA Remarques en vrac:
	18	cAA--------------------
	19	cAA 1/ le call phytrac se fait avec nqmax-2 donc nous avons bien
	20	cAA les vrais traceurs (nbtr) dans phytrac (pas la vapeur ni eau liquide)
	21	cAA 2/ Le choix du radon et du pb se fait juste avec un data
	22	cAA (peu propre). Peut-etre pourrait-on prevoir dans l'avenir
	23	cAA une variable "type de traceur"
	24	c======================================================================
	25	#include "YOMCST.h"
	26	#include "dimensions.h"
	27	#include "dimphy.h"
	28	#include "indicesol.h"
	29	#include "control.h"
[53]	30	#include "temps.h"
[2]	31	c======================================================================
	32
	33	c Arguments:
	34	c
	35	c EN ENTREE:
	36	c ==========
	37	c
	38	c divers:
	39	c -------
	40	c
	41	integer nlon ! nombre de points horizontaux
	42	integer nlev ! nombre de couches verticales
	43	integer nqmax ! nombre de traceurs auxquels on applique la physique
	44	real pdtphys ! pas d'integration pour la physique (seconde)
	45	real t_seri(nlon,nlev) ! temperature
	46	real tr_seri(nlon,nlev,nbtr) ! traceur
	47	real paprs(nlon,nlev+1) ! pression pour chaque inter-couche (en Pa)
	48	real pplay(nlon,nlev) ! pression pour le mileu de chaque couche (en Pa)
	49	real presnivs(klev) ! pressions approximat. des milieux couches ( en PA)
	50	real znivsig(klev) ! niveaux sigma
	51	real paire(klon)
	52	real pphis(klon)
	53	logical debutphy ! le flag de l'initialisation de la physique
	54	integer ll
	55	c
	56	cAA Rem : nbtr : nombre de vrais traceurs est defini dans dimphy.h
	57	c
	58	c convection:
	59	c -----------
	60	c
	61	REAL pmfu(nlon,nlev) ! flux de masse dans le panache montant
	62	REAL pmfd(nlon,nlev) ! flux de masse dans le panache descendant
	63	REAL pen_u(nlon,nlev) ! flux entraine dans le panache montant
	64	REAL pde_u(nlon,nlev) ! flux detraine dans le panache montant
	65	REAL pen_d(nlon,nlev) ! flux entraine dans le panache descendant
	66	REAL pde_d(nlon,nlev) ! flux detraine dans le panache descendant
	67	c
	68	c Couche limite:
	69	c --------------
	70	c
	71	REAL coefh(nlon,nlev) ! coeff melange CL
	72	REAL yu1(nlon) ! vents au premier niveau
	73	REAL yv1(nlon) ! vents au premier niveau
	74	REAL xlat(nlon) ! latitudes pour chaque point
	75	REAL xlon(nlon) ! longitudes pour chaque point
	76
	77	c
	78	c Lessivage:
	79	c ----------
	80	c
	81	c pour le ON-LINE
	82	c
	83	REAL frac_impa(nlon,nlev) ! fraction d'aerosols impactes
	84	REAL frac_nucl(nlon,nlev) ! fraction d'aerosols nuclees
	85	c
	86	cAA
	87	cAA Arguments necessaires pour les sources et puits de traceur:
	88	cAA ----------------
	89	cAA
	90	real ftsol(nlon,nbsrf) ! Temperature du sol (surf)(Kelvin)
	91	real pctsrf(nlon,nbsrf) ! Pourcentage de sol f(nature du sol)
	92
	93	cAA ----------------------------
	94	cAA VARIABLES LOCALES TRACEURS
	95	cAA ----------------------------
	96	cAA
	97	cAA Sources et puits des traceurs:
	98	cAA ------------------------------
	99	cAA
	100	cAA Pour l'instant seuls les cas du rn et du pb ont ete envisages.
	101
	102	REAL source(klon) ! a voir lorsque le flux est prescrit
	103	REAL topuits(klev,nbtr) ! a voir
	104	cAA
	105	cAA Pour la source de radon et son reservoir de sol
	106	cAA ................................................
	107
	108	REAL trs(klon,nbtr) ! Conc. radon ds le sol
	109	SAVE trs
	110
	111	REAL masktr(klon,nbtr) ! Masque reservoir de sol traceur
	112	c Masque de l'echange avec la surface
	113	c (1 = reservoir) ou (possible => 1 )
	114	SAVE masktr
	115	REAL fshtr(klon,nbtr) ! Flux surfacique dans le reservoir de sol
	116	SAVE fshtr
	117	REAL hsoltr(nbtr) ! Epaisseur equivalente du reservoir de sol
	118	SAVE hsoltr
	119	REAL tautr(nbtr) ! Constante de decroissance radioactive
	120	SAVE tautr
	121	REAL vdeptr(nbtr) ! Vitesse de depot sec dans la couche Brownienne
	122	SAVE vdeptr
	123	REAL scavtr(nbtr) ! Coefficient de lessivage
	124	SAVE scavtr
	125	cAA
	126	CHARACTER*1 itn
	127	C maf ioipsl
	128	CHARACTER*2 str2
	129	INTEGER nhori, nvert
	130	REAL zsto, zout, zjulian
	131	INTEGER nid_tra
	132	SAVE nid_tra
[53]	133	INTEGER ndex2d(iim(jjm+1)),ndex3d(iim(jjm+1)*klev)
[2]	134	REAL zx_tmp_2d(iim,jjm+1), zx_tmp_3d(iim,jjm+1,klev)
	135	REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1)
	136	INTEGER itra
	137	SAVE itra ! compteur pour la physique
	138	C
	139	C Variables liees a l'ecriture de la bande histoire : phytrac.nc
	140	c
	141	INTEGER ecrit_tra
	142	SAVE ecrit_tra
[53]	143	logical ok_sync
	144	parameter (ok_sync = .true.)
[2]	145	C
	146	C nature du traceur
	147	c
	148	logical aerosol(nbtr) ! Nature du traceur
	149	c ! aerosol(it) = true => aerosol
	150	c ! aerosol(it) = false => gaz
	151	c ! nat_trac(it) = 1. aerosolc
	152	logical clsol(nbtr) ! clsol(it) = true => CL sol calculee
	153	logical radio(nbtr) ! radio(it)=true => decroisssance radioactive
	154	save aerosol,clsol,radio
	155	C
	156	c======================================================================
	157	c
	158	c Declaration des procedures appelees
	159	c
	160	c--modif convection tiedtke
	161	INTEGER i, k, it
	162
	163	REAL delp(klon,klev)
	164	c--end modif
	165	c
	166	c Variables liees a l'ecriture de la bande histoire physique
	167	c
	168	c Variables locales pour effectuer les appels en serie
	169	c----------------------------------------------------
	170	c
	171	REAL d_tr(klon,klev), d_trs(klon) ! tendances de traceurs
	172	REAL d_tr_cl(klon,klev) ! tendance de traceurs couche limite
	173	REAL d_tr_cv(klon,klev) ! tendance de traceurs convection
	174	REAL d_tr_dec(klon,klev,nbtr) ! la tendance de la decroissance
	175	c ! radioactive du rn - > pb
	176	REAL d_tr_lessi_impa(klon,klev,nbtr) ! la tendance du lessivage
	177	c ! par impaction
	178	REAL d_tr_lessi_nucl(klon,klev,nbtr) ! la tendance du lessivage
	179	c ! par nucleation
	180	REAL fluxrn(klon,klev)
	181	REAL fluxpb(klon,klev)
	182	REAL pbimpa(klon,klev)
	183	REAL pbnucl(klon,klev)
	184	REAL rn(klon,klev)
	185	REAL pb(klon,klev)
	186	REAL flestottr(klon,klev,nbtr) ! flux de lessivage
	187	c ! dans chaque couche
	188
	189	C
	190	character*20 modname
	191	character*80 abort_message
	192	c
	193	c Controles
	194	c-------------
	195	logical first,couchelimite,convection,lessivage,sorties,
	196	s rnpb,inirnpb
	197	save first,couchelimite,convection,lessivage,sorties,
	198	s inirnpb
	199	data first,couchelimite,convection,lessivage,sorties
	200	s /.true.,.true.,.true.,.true.,.true./
	201	c
	202	c======================================================================
	203	c ecrit_tra = NINT(86400./pdtphys *1.0) ! tous les jours
	204	modname='phytrac'
	205
[53]	206	c print*,'DANS PHYTRAC debutphy=',debutphy
[2]	207
	208	ecrit_tra = NINT(86400./pdtphys *ecritphy)
[53]	209	zsto = pdtphys
	210	zout = pdtphys * FLOAT(ecrit_tra)
[2]	211	if (debutphy) then
	212
	213	if(nbtr.lt.nqmax) then
	214	print*,'NQMAX=',nqmax
	215	print*,'NBTR=',nbtr
	216	abort_message='See above'
	217	call abort_gcm(modname,abort_message,1)
	218	endif
	219
	220	inirnpb=rnpb
	221	PRINT*, 'La frequence de sortie traceurs est ', ecrit_tra
	222	itra=0
	223	C
[53]	224	CALL ymds2ju(anne_ini, 1, 1, 0.0, zjulian)
	225	zjulian = zjulian + day_ini
[2]	226	c
	227	CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlon,zx_lon)
	228	DO i = 1, iim
	229	zx_lon(i,1) = xlon(i+1)
	230	zx_lon(i,jjm+1) = xlon(i+1)
	231	ENDDO
	232	DO ll=1,klev
	233	znivsig(ll)=float(ll)
	234	ENDDO
	235	CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlat,zx_lat)
	236	CALL histbeg("histrac", iim,zx_lon, jjm+1,zx_lat,
	237	. 1,iim,1,jjm+1, 0, zjulian, pdtphys,
	238	. nhori, nid_tra)
	239	call histvert(nid_tra, 'sig_s', 'Niveaux sigma','-',
	240	. klev, znivsig, nvert)
	241	C
	242	C CALL histvert(nid_tra, "presnivs", "Vertical levels", "mb",
	243	C . klev, presnivs, nvert)
	244	c
	245	CALL histdef(nid_tra, "phis", "Surface geop. height", "-",
	246	. iim,jjm+1,nhori, 1,1,1, -99, 32,
	247	. "once", zsto,zout)
	248	c
	249	CALL histdef(nid_tra, "aire", "Grid area", "-",
	250	. iim,jjm+1,nhori, 1,1,1, -99, 32,
	251	. "once", zsto,zout)
	252	c
	253	DO it=1,nqmax
	254	IF (it.LE.99) THEN
	255	WRITE(str2,'(i2.2)') it
	256	C champ 3D
	257	CALL histdef(nid_tra, "tr"//str2, "Tracer No."//str2, "U/kga",
	258	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	259	. "ave(X)", zsto,zout)
	260	IF (lessivage) THEN
	261	CALL histdef(nid_tra, "fl"//str2, "Tracer Flux No."//str2,
	262	. "U/m2/s",iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	263	. "ave(X)", zsto,zout)
	264	ENDIF
	265	ELSE
	266	PRINT*, "Trop de traceurs"
	267	CALL abort
	268	ENDIF
	269	ENDDO
	270	CALL histend(nid_tra)
[53]	271
[2]	272	c======================================================================
	273	c Initialisation de certaines variables pour le Rn et le Pb
	274	c======================================================================
	275
	276	c Initialisation du traceur dans le sol (couche limite radonique)
	277	c
[53]	278	c print*,'valeur de debut dans phytrac :',debutphy
[2]	279	do it=1,nqmax
	280	do i=1,klon
	281	trs(i,it) = 0.
	282	enddo
	283	END DO
	284	c Initialisation de la fraction d'aerosols lessivee
	285	c
	286	DO it=1,nqmax
	287	DO k = 1, nlev
	288	DO i = 1, klon
	289	d_tr_lessi_impa(i,k,it) = 0.0
	290	d_tr_lessi_nucl(i,k,it) = 0.0
	291	ENDDO
	292	ENDDO
	293	ENDDO
	294	c
	295	c Initialisation de la nature des traceurs
	296	c
	297	DO it = 1, nqmax
	298	aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut
	299	radio(it) = .FALSE. ! Par defaut pas de passage par radiornpb
	300	clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit
	301	ENDDO
	302	c
	303	ENDIF ! fin debutphy
	304	c Initialisation du traceur dans le sol (couche limite radonique)
	305	if(inirnpb) THEN
	306	c
	307	radio(1)= .true.
	308	radio(2)= .true.
	309	clsol(1)= .true.
	310	clsol(2)= .true.
	311	aerosol(2) = .TRUE. ! le Pb est un aerosol
	312	c
	313	call initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr
	314	. ,vdeptr,scavtr)
	315	inirnpb=.false.
	316	endif
	317	c======================================================================
	318	c Calcul de l'effet de la convection
	319	c======================================================================
	320
	321	if (convection) then
	322
[53]	323	c print*,'Pas de temps dans phytrac : ',pdtphys
[2]	324	DO it=1, nqmax
	325	CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
	326	. pplay, paprs, tr_seri(1,1,it), d_tr_cv)
	327	DO k = 1, nlev
	328	DO i = 1, klon
	329	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k)
	330	ENDDO
	331	ENDDO
	332	WRITE(itn,'(i1)') it
	333	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it='//itn)
	334	ENDDO
[53]	335	c print*,'apres nflxtr'
[2]	336
	337
	338	endif ! convection
	339
	340	c======================================================================
	341	c Calcul de l'effet de la couche limite
	342	c======================================================================
	343
[53]	344	c print*,'avant couchelimite'
[2]	345	if (couchelimite) then
	346
	347	DO k = 1, nlev
	348	DO i = 1, klon
	349	delp(i,k) = paprs(i,k)-paprs(i,k+1)
	350	ENDDO
	351	ENDDO
	352
	353	C maf modif pour tenir compte du cas rnpb + traceur
	354	DO it=1, nqmax
[53]	355	c print *,'it',it,clsol(it)
[2]	356	if (clsol(it)) then ! couche limite avec quantite dans le sol calculee
	357	CALL cltracrn(it, pdtphys, yu1, yv1,
	358	e coefh,t_seri,ftsol,pctsrf,
	359	e tr_seri(1,1,it),trs(1,it),
	360	e paprs, pplay, delp,
	361	e masktr(1,it),fshtr(1,it),hsoltr(it),
	362	e tautr(it),vdeptr(it),
	363	e xlat,
	364	s d_tr_cl,d_trs)
	365	DO k = 1, nlev
	366	DO i = 1, klon
	367	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k)
	368	ENDDO
	369	ENDDO
	370	c
	371	c Traceur ds sol
	372	c
	373	DO i = 1, klon
	374	trs(i,it) = trs(i,it) + d_trs(i)
	375	END DO
	376	C
	377	C maf provisoire suppression des prints
	378	C WRITE(itn,'(i1)') it
	379	C CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
	380	else ! couche limite avec flux prescrit
	381	Cmaf provisoire source / traceur a creer
	382	DO i=1, klon
	383	source(i) = 0.0 ! pas de source, pour l'instant
	384	ENDDO
	385	C
	386	CALL cltrac(pdtphys, coefh,t_seri,
	387	s tr_seri(1,1,it), source,
	388	e paprs, pplay, delp,
	389	s d_tr )
	390	DO k = 1, nlev
	391	DO i = 1, klon
	392	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr(i,k)
	393	ENDDO
	394	ENDDO
	395	Cmaf provisoire suppression des prints
	396	Cmaf WRITE(itn,'(i1)') it
	397	cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
	398	endif
	399	ENDDO
	400	c
	401	endif ! couche limite
	402
[53]	403	c print*,'apres couchelimite'
[2]	404
	405	c======================================================================
	406	c Calcul de l'effet du puits radioactif
	407	c======================================================================
	408
	409	C MAF il faudrait faire une modification pour passer dans radiornpb
	410	c si radio=true mais pour l'instant radiornpb propre au cas rnpb
	411	if(rnpb) then
	412	call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
	413	C
	414	DO it=1,nqmax
	415	if(radio(it)) then
	416	DO k = 1, nlev
	417	DO i = 1, klon
	418	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
	419	ENDDO
	420	ENDDO
	421	WRITE(itn,'(i1)') it
	422	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn)
	423	endif
	424	ENDDO
	425	c
	426	endif ! rnpb decroissance radioactive
	427	C
	428	c======================================================================
	429	c Calcul de l'effet de la precipitation
	430	c======================================================================
	431
	432	IF (lessivage) THEN
	433
[53]	434	c print*,'avant lessivage'
[2]	435
	436	DO it = 1, nqmax
	437	DO k = 1, nlev
	438	DO i = 1, klon
	439	d_tr_lessi_nucl(i,k,it) = 0.0
	440	d_tr_lessi_impa(i,k,it) = 0.0
	441	flestottr(i,k,it) = 0.0
	442	ENDDO
	443	ENDDO
	444	ENDDO
	445	c
	446	c tendance des aerosols nuclees et impactes
	447	c
	448	DO it = 1, nqmax
	449	IF (aerosol(it)) THEN
	450	DO k = 1, nlev
	451	DO i = 1, klon
	452	d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +
	453	s ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
	454	d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +
	455	s ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
	456	ENDDO
	457	ENDDO
	458	ENDIF
	459	ENDDO
	460	c
	461	c Mises a jour des traceurs + calcul des flux de lessivage
	462	c Mise a jour due a l'impaction et a la nucleation
	463	c
	464	DO it = 1, nqmax
	465	IF (aerosol(it)) THEN
	466	DO k = 1, nlev
	467	DO i = 1, klon
	468	tr_seri(i,k,it) = tr_seri(i,k,it) *
	469	s ( frac_impa(i,k) + frac_nucl(i,k) - 1. )
	470	ENDDO
	471	ENDDO
	472	ENDIF
	473	ENDDO
	474	c
	475	c Flux lessivage total
	476	c
	477	DO it = 1, nqmax
	478	DO k = 1, nlev
	479	DO i = 1, klon
	480	flestottr(i,k,it) = flestottr(i,k,it) -
	481	s ( d_tr_lessi_nucl(i,k,it) +
	482	s d_tr_lessi_impa(i,k,it) ) *
	483	s ( paprs(i,k)-paprs(i,k+1) ) /
	484	s (RG * pdtphys)
	485	ENDDO
	486	ENDDO
	487	c
	488	Cmaf suppression provisoire
	489	Cmaf WRITE(itn,'(i1)') it
	490	Cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
	491	ENDDO
	492	c
[53]	493	c print*,'apres lessivage'
[2]	494	ENDIF
	495	Cc
	496	DO k = 1, nlev
	497	DO i = 1, klon
	498	fluxrn(i,k) = flestottr(i,k,1)
	499	fluxpb(i,k) = flestottr(i,k,2)
	500	rn(i,k) = tr_seri(i,k,1)
	501	pb(i,k) = tr_seri(i,k,2)
	502	pbnucl(i,k)=d_tr_lessi_nucl(i,k,2)
	503	pbimpa(i,k)=d_tr_lessi_impa(i,k,2)
	504	ENDDO
	505	ENDDO
	506	itra=itra+1
[53]	507
	508	C
	509	C Sorties IOIPSL
	510	ndex2d = 0
	511	ndex3d = 0
	512	c
	513	c write(,)'sorties ioipsl phytrac',zsto,zout
	514	CALL gr_fi_ecrit(1,klon,iim,jjm+1,pphis,zx_tmp_2d)
	515	CALL histwrite(nid_tra,"phis",itra,zx_tmp_2d,iim*(jjm+1),ndex2d)
	516	C
	517	CALL gr_fi_ecrit(1,klon,iim,jjm+1,paire,zx_tmp_2d)
	518	CALL histwrite(nid_tra,"aire",itra,zx_tmp_2d,iim*(jjm+1),ndex2d)
[2]	519	DO it=1,nqmax
	520	IF (it.LE.99) THEN
	521	WRITE(str2,'(i2.2)') it
	522	C champs 3D
	523	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,tr_seri(1,1,it),zx_tmp_3d)
	524	CALL histwrite(nid_tra,"tr"//str2,itra,zx_tmp_3d,
[53]	525	. iim(jjm+1)klev,ndex3d)
[2]	526	IF (lessivage) THEN
	527	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,flestottr(1,1,it),zx_tmp_3d)
	528	CALL histwrite(nid_tra,"fl"//str2,itra,zx_tmp_3d,
[53]	529	. iim(jjm+1)klev,ndex3d)
[2]	530	ENDIF
	531	ELSE
	532	PRINT*, "Trop de traceurs"
	533	CALL abort
	534	ENDIF
	535	ENDDO
[53]	536	if (ok_sync) call histsync(nid_tra)
[2]	537
	538	RETURN
	539	END

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

Download in other formats:

Original Format