Context Navigation

phytrac.F @ 199

Last change on this file since 199 was 199, checked in by lmdz, 24 years ago
Introduction Kerry Emmanuel ratqs interactif mise a jour offline FH LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 18.1 KB

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

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

Download in other formats:

Original Format