Context Navigation

phytrac.F @ 332

Last change on this file since 332 was 209, checked in by lmdz, 24 years ago
Manque deux include JQuaas LF
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 26.0 KB

Rev	Line
[199]	1	c
[207]	2	c $Header$
[199]	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"
[209]	32	#include "temps.h"
	33	#include "control.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 paire(klon)
	54	real pphis(klon)
	55	logical debutphy ! le flag de l'initialisation de la physique
[199]	56	logical lafin ! le flag de la fin de la physique
	57
[2]	58	integer ll
	59	c
	60	cAA Rem : nbtr : nombre de vrais traceurs est defini dans dimphy.h
	61	c
	62	c convection:
	63	c -----------
	64	c
	65	REAL pmfu(nlon,nlev) ! flux de masse dans le panache montant
	66	REAL pmfd(nlon,nlev) ! flux de masse dans le panache descendant
	67	REAL pen_u(nlon,nlev) ! flux entraine dans le panache montant
	68	REAL pde_u(nlon,nlev) ! flux detraine dans le panache montant
	69	REAL pen_d(nlon,nlev) ! flux entraine dans le panache descendant
	70	REAL pde_d(nlon,nlev) ! flux detraine dans le panache descendant
	71	c
	72	c Couche limite:
	73	c --------------
	74	c
	75	REAL coefh(nlon,nlev) ! coeff melange CL
	76	REAL yu1(nlon) ! vents au premier niveau
	77	REAL yv1(nlon) ! vents au premier niveau
	78	REAL xlat(nlon) ! latitudes pour chaque point
	79	REAL xlon(nlon) ! longitudes pour chaque point
	80
	81	c
	82	c Lessivage:
	83	c ----------
	84	c
	85	c pour le ON-LINE
	86	c
	87	REAL frac_impa(nlon,nlev) ! fraction d'aerosols impactes
	88	REAL frac_nucl(nlon,nlev) ! fraction d'aerosols nuclees
	89	c
	90	cAA
	91	cAA Arguments necessaires pour les sources et puits de traceur:
	92	cAA ----------------
	93	cAA
	94	real ftsol(nlon,nbsrf) ! Temperature du sol (surf)(Kelvin)
	95	real pctsrf(nlon,nbsrf) ! Pourcentage de sol f(nature du sol)
[204]	96	c abder
	97	real pftsol1(nlon),pftsol2(nlon),pftsol3(nlon),pftsol4(nlon)
	98	real ppsrf1(nlon),ppsrf2(nlon),ppsrf3(nlon),ppsrf4(nlon)
	99	c fin
[2]	100	cAA ----------------------------
	101	cAA VARIABLES LOCALES TRACEURS
	102	cAA ----------------------------
	103	cAA
	104	cAA Sources et puits des traceurs:
	105	cAA ------------------------------
	106	cAA
	107	cAA Pour l'instant seuls les cas du rn et du pb ont ete envisages.
	108
	109	REAL source(klon) ! a voir lorsque le flux est prescrit
	110	REAL topuits(klev,nbtr) ! a voir
	111	cAA
	112	cAA Pour la source de radon et son reservoir de sol
	113	cAA ................................................
	114
	115	REAL trs(klon,nbtr) ! Conc. radon ds le sol
	116	SAVE trs
	117
	118	REAL masktr(klon,nbtr) ! Masque reservoir de sol traceur
	119	c Masque de l'echange avec la surface
	120	c (1 = reservoir) ou (possible => 1 )
	121	SAVE masktr
	122	REAL fshtr(klon,nbtr) ! Flux surfacique dans le reservoir de sol
	123	SAVE fshtr
	124	REAL hsoltr(nbtr) ! Epaisseur equivalente du reservoir de sol
	125	SAVE hsoltr
	126	REAL tautr(nbtr) ! Constante de decroissance radioactive
	127	SAVE tautr
	128	REAL vdeptr(nbtr) ! Vitesse de depot sec dans la couche Brownienne
	129	SAVE vdeptr
	130	REAL scavtr(nbtr) ! Coefficient de lessivage
	131	SAVE scavtr
	132	cAA
	133	CHARACTER*1 itn
	134	C maf ioipsl
	135	CHARACTER*2 str2
	136	INTEGER nhori, nvert
	137	REAL zsto, zout, zjulian
	138	INTEGER nid_tra
	139	SAVE nid_tra
[204]	140	c REAL x(klon,klev,nbtr+2) ! traceurs
	141	INTEGER ndex(1)
[2]	142	REAL zx_tmp_2d(iim,jjm+1), zx_tmp_3d(iim,jjm+1,klev)
	143	REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1)
	144	INTEGER itra
	145	SAVE itra ! compteur pour la physique
	146	C
	147	C Variables liees a l'ecriture de la bande histoire : phytrac.nc
	148	c
	149	INTEGER ecrit_tra
	150	SAVE ecrit_tra
[53]	151	logical ok_sync
	152	parameter (ok_sync = .true.)
[2]	153	C
	154	C nature du traceur
	155	c
	156	logical aerosol(nbtr) ! Nature du traceur
	157	c ! aerosol(it) = true => aerosol
	158	c ! aerosol(it) = false => gaz
	159	c ! nat_trac(it) = 1. aerosolc
	160	logical clsol(nbtr) ! clsol(it) = true => CL sol calculee
	161	logical radio(nbtr) ! radio(it)=true => decroisssance radioactive
	162	save aerosol,clsol,radio
	163	C
	164	c======================================================================
	165	c
	166	c Declaration des procedures appelees
	167	c
	168	c--modif convection tiedtke
[204]	169	INTEGER i, k, it,itap
	170	save itap
[2]	171	REAL delp(klon,klev)
	172	c--end modif
	173	c
	174	c Variables liees a l'ecriture de la bande histoire physique
	175	c
	176	c Variables locales pour effectuer les appels en serie
	177	c----------------------------------------------------
	178	c
	179	REAL d_tr(klon,klev), d_trs(klon) ! tendances de traceurs
	180	REAL d_tr_cl(klon,klev) ! tendance de traceurs couche limite
	181	REAL d_tr_cv(klon,klev) ! tendance de traceurs convection
	182	REAL d_tr_dec(klon,klev,nbtr) ! la tendance de la decroissance
	183	c ! radioactive du rn - > pb
	184	REAL d_tr_lessi_impa(klon,klev,nbtr) ! la tendance du lessivage
	185	c ! par impaction
	186	REAL d_tr_lessi_nucl(klon,klev,nbtr) ! la tendance du lessivage
	187	c ! par nucleation
	188	REAL fluxrn(klon,klev)
	189	REAL fluxpb(klon,klev)
	190	REAL pbimpa(klon,klev)
	191	REAL pbnucl(klon,klev)
	192	REAL rn(klon,klev)
	193	REAL pb(klon,klev)
	194	REAL flestottr(klon,klev,nbtr) ! flux de lessivage
	195	c ! dans chaque couche
	196
	197	C
	198	character*20 modname
	199	character*80 abort_message
	200	c
	201	c Controles
	202	c-------------
	203	logical first,couchelimite,convection,lessivage,sorties,
	204	s rnpb,inirnpb
	205	save first,couchelimite,convection,lessivage,sorties,
	206	s inirnpb
	207	data first,couchelimite,convection,lessivage,sorties
	208	s /.true.,.true.,.true.,.true.,.true./
	209	c
	210	c======================================================================
	211	c ecrit_tra = NINT(86400./pdtphys *1.0) ! tous les jours
	212	modname='phytrac'
	213
[53]	214	c print*,'DANS PHYTRAC debutphy=',debutphy
[2]	215
	216	if (debutphy) then
	217
[204]	218	print*,'dans phytrac ',pdtphys,ecritphy,ecrit_tra
	219	ecrit_tra = NINT(86400./pdtphys/2.) ! tous les 12H
	220	c ecrit_tra = NINT(86400./pdtphys) ! tous les 24H
	221
[2]	222	if(nbtr.lt.nqmax) then
	223	print*,'NQMAX=',nqmax
	224	print*,'NBTR=',nbtr
	225	abort_message='See above'
	226	call abort_gcm(modname,abort_message,1)
	227	endif
	228
	229	inirnpb=rnpb
	230	PRINT*, 'La frequence de sortie traceurs est ', ecrit_tra
	231	itra=0
[204]	232	itap=0
[2]	233	C
[53]	234	CALL ymds2ju(anne_ini, 1, 1, 0.0, zjulian)
	235	zjulian = zjulian + day_ini
[2]	236	c
	237	CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlon,zx_lon)
	238	DO i = 1, iim
	239	zx_lon(i,1) = xlon(i+1)
	240	zx_lon(i,jjm+1) = xlon(i+1)
	241	ENDDO
[209]	242	c DO ll=1,klev
	243	c znivsig(ll)=float(ll)
	244	c ENDDO
[2]	245	CALL gr_fi_ecrit(1,klon,iim,jjm+1,xlat,zx_lat)
	246	CALL histbeg("histrac", iim,zx_lon, jjm+1,zx_lat,
	247	. 1,iim,1,jjm+1, 0, zjulian, pdtphys,
	248	. nhori, nid_tra)
[204]	249	CALL histvert(nid_tra, "presnivs", "Vertical levels", "mb",
	250	. klev, presnivs, nvert)
	251	zsto = pdtphys
	252	zout = pdtphys * FLOAT(ecrit_tra)
[2]	253	c
	254	CALL histdef(nid_tra, "phis", "Surface geop. height", "-",
	255	. iim,jjm+1,nhori, 1,1,1, -99, 32,
	256	. "once", zsto,zout)
	257	c
	258	CALL histdef(nid_tra, "aire", "Grid area", "-",
	259	. iim,jjm+1,nhori, 1,1,1, -99, 32,
	260	. "once", zsto,zout)
[204]	261
	262	goto 666
	263	CALL histdef(nid_tra, "pyu1", "Vent niv 1", "-",
	264	. iim,jjm+1,nhori, 1,1,1, -99, 32,
	265	. "inst(X)", zsto,zout)
	266
	267	CALL histdef(nid_tra, "pyv1", "Vent niv 1", "-",
	268	. iim,jjm+1,nhori, 1,1,1, -99, 32,
	269	. "inst(X)", zsto,zout)
	270	CALL histdef(nid_tra, "psrf1", "nature sol", "-",
	271	. iim,jjm+1,nhori, 1,1,1, -99, 32,
	272	. "inst(X)", zsto,zout)
	273	CALL histdef(nid_tra, "psrf2", "nature sol", "-",
	274	. iim,jjm+1,nhori, 1,1,1, -99, 32,
	275	. "inst(X)", zsto,zout)
	276	CALL histdef(nid_tra, "psrf3", "nature sol", "-",
	277	. iim,jjm+1,nhori, 1,1,1, -99, 32,
	278	. "inst(X)", zsto,zout)
	279	CALL histdef(nid_tra, "psrf4", "nature sol", "-",
	280	. iim,jjm+1,nhori, 1,1,1, -99, 32,
	281	. "inst(X)", zsto,zout)
	282	CALL histdef(nid_tra, "ftsol1", "temper sol", "-",
	283	. iim,jjm+1,nhori, 1,1,1, -99, 32,
	284	. "inst(X)", zsto,zout)
	285	CALL histdef(nid_tra, "ftsol2", "temper sol", "-",
	286	. iim,jjm+1,nhori, 1,1,1, -99, 32,
	287	. "inst(X)", zsto,zout)
	288	CALL histdef(nid_tra, "ftsol3", "temper sol", "-",
	289	. iim,jjm+1,nhori, 1,1,1, -99, 32,
	290	. "inst", zsto,zout)
	291	CALL histdef(nid_tra, "ftsol4", "temper sol", "-",
	292	. iim,jjm+1,nhori, 1,1,1, -99, 32,
	293	. "inst(X)", zsto,zout)
	294	CALL histdef(nid_tra, "pplay", "flux u mont","-",
	295	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	296	. "inst(X)", zsto,zout)
	297	CALL histdef(nid_tra, "t", "flux u mont","-",
	298	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	299	. "inst(X)", zsto,zout)
	300	CALL histdef(nid_tra, "mfu", "flux u mont","-",
	301	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	302	. "ave(X)", zsto,zout)
	303	CALL histdef(nid_tra, "mfd", "flux u decen","-",
	304	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	305	. "ave(X)", zsto,zout)
	306	CALL histdef(nid_tra, "en_u", "flux u mont","-",
	307	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	308	. "ave(X)", zsto,zout)
	309	CALL histdef(nid_tra, "en_d", "flux u mont","-",
	310	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	311	. "ave(X)", zsto,zout)
	312	CALL histdef(nid_tra, "de_u", "flux u mont","-",
	313	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	314	. "ave(X)", zsto,zout)
	315	CALL histdef(nid_tra, "de_d", "flux u mont","-",
	316	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	317	. "ave(X)", zsto,zout)
	318	CALL histdef(nid_tra, "coefh", "turbulent coef","-",
	319	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	320	. "ave(X)", zsto,zout)
	321
	322	666 continue
[2]	323	c
	324	DO it=1,nqmax
	325	IF (it.LE.99) THEN
	326	WRITE(str2,'(i2.2)') it
	327	C champ 3D
	328	CALL histdef(nid_tra, "tr"//str2, "Tracer No."//str2, "U/kga",
	329	. iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	330	. "ave(X)", zsto,zout)
	331	IF (lessivage) THEN
	332	CALL histdef(nid_tra, "fl"//str2, "Tracer Flux No."//str2,
	333	. "U/m2/s",iim,jjm+1,nhori, klev,1,klev,nvert, 32,
	334	. "ave(X)", zsto,zout)
	335	ENDIF
	336	ELSE
	337	PRINT*, "Trop de traceurs"
	338	CALL abort
	339	ENDIF
	340	ENDDO
	341	CALL histend(nid_tra)
[204]	342	ndex(1) = 0
	343	c
	344	i = NINT(zout/zsto)
	345	CALL gr_fi_ecrit(1,klon,iim,jjm+1,pphis,zx_tmp_2d)
	346	CALL histwrite(nid_tra,"phis",i,zx_tmp_2d,iim*(jjm+1),ndex)
	347	C
	348	i = NINT(zout/zsto)
	349	CALL gr_fi_ecrit(1,klon,iim,jjm+1,paire,zx_tmp_2d)
	350	CALL histwrite(nid_tra,"aire",i,zx_tmp_2d,iim*(jjm+1),ndex)
[53]	351
[2]	352	c======================================================================
	353	c Initialisation de certaines variables pour le Rn et le Pb
	354	c======================================================================
	355
	356	c Initialisation du traceur dans le sol (couche limite radonique)
	357	c
[53]	358	c print*,'valeur de debut dans phytrac :',debutphy
[2]	359	do it=1,nqmax
	360	do i=1,klon
	361	trs(i,it) = 0.
	362	enddo
	363	END DO
[199]	364
	365	open (99,file='starttrac',status='old',
	366	. err=999,form='formatted')
	367	read(99,*) (trs(i,1),i=1,klon)
	368	999 close(99)
	369	print*, 'apres starttrac'
	370
[2]	371	c Initialisation de la fraction d'aerosols lessivee
	372	c
	373	DO it=1,nqmax
	374	DO k = 1, nlev
	375	DO i = 1, klon
	376	d_tr_lessi_impa(i,k,it) = 0.0
	377	d_tr_lessi_nucl(i,k,it) = 0.0
	378	ENDDO
	379	ENDDO
	380	ENDDO
	381	c
	382	c Initialisation de la nature des traceurs
	383	c
	384	DO it = 1, nqmax
	385	aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut
	386	radio(it) = .FALSE. ! Par defaut pas de passage par radiornpb
	387	clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit
	388	ENDDO
	389	c
	390	ENDIF ! fin debutphy
	391	c Initialisation du traceur dans le sol (couche limite radonique)
	392	if(inirnpb) THEN
	393	c
	394	radio(1)= .true.
	395	radio(2)= .true.
	396	clsol(1)= .true.
	397	clsol(2)= .true.
	398	aerosol(2) = .TRUE. ! le Pb est un aerosol
	399	c
	400	call initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr
	401	. ,vdeptr,scavtr)
	402	inirnpb=.false.
	403	endif
[204]	404	if(nqmax.gt.2) aerosol(3)=.true.
	405
	406
	407	c abder
	408	goto 777
	409	do i=1,nlon
	410	pftsol1(i) = ftsol(i,1)
	411	pftsol2(i) = ftsol(i,2)
	412	pftsol3(i) = ftsol(i,3)
	413	pftsol4(i) = ftsol(i,4)
	414
	415	ppsrf1(i) = pctsrf(i,1)
	416	ppsrf2(i) = pctsrf(i,2)
	417	ppsrf3(i) = pctsrf(i,3)
	418	ppsrf4(i) = pctsrf(i,4)
	419
	420	enddo
	421	ndex(1)=0
	422	itap=itap+1
	423	CALL gr_fi_ecrit(1,klon,iim,jjm+1,yu1,zx_tmp_2d)
	424	CALL histwrite(nid_tra,"pyu1",itap,zx_tmp_2d,
	425	s iim*(jjm+1),ndex)
	426
	427	CALL gr_fi_ecrit(1,klon,iim,jjm+1,yv1,zx_tmp_2d)
	428	CALL histwrite(nid_tra,"pyv1",itap,zx_tmp_2d,
	429	s iim*(jjm+1),ndex)
	430
	431	CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol1,zx_tmp_2d)
	432	CALL histwrite(nid_tra,"ftsol1",itap,zx_tmp_2d,
	433	s iim*(jjm+1),ndex)
	434
	435	CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol2,zx_tmp_2d)
	436	CALL histwrite(nid_tra,"ftsol2",itap,zx_tmp_2d,
	437	s iim*(jjm+1),ndex)
	438
	439	CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol3,zx_tmp_2d)
	440	CALL histwrite(nid_tra,"ftsol3",itap,zx_tmp_2d,
	441	s iim*(jjm+1),ndex)
	442
	443	CALL gr_fi_ecrit(1,klon,iim,jjm+1,pftsol4,zx_tmp_2d)
	444	CALL histwrite(nid_tra,"ftsol4",itap,zx_tmp_2d,
	445	s iim*(jjm+1),ndex)
	446
	447	CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf1,zx_tmp_2d)
	448	CALL histwrite(nid_tra,"psrf1",itap,zx_tmp_2d,
	449	s iim*(jjm+1),ndex)
	450
	451	CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf2,zx_tmp_2d)
	452	CALL histwrite(nid_tra,"psrf2",itap,zx_tmp_2d,
	453	s iim*(jjm+1),ndex)
	454
	455	CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf3,zx_tmp_2d)
	456	CALL histwrite(nid_tra,"psrf3",itap,zx_tmp_2d,
	457	s iim*(jjm+1),ndex)
	458
	459	CALL gr_fi_ecrit(1,klon,iim,jjm+1,ppsrf4,zx_tmp_2d)
	460	CALL histwrite(nid_tra,"psrf4",itap,zx_tmp_2d,
	461	s iim*(jjm+1),ndex)
	462	777 continue
[2]	463	c======================================================================
	464	c Calcul de l'effet de la convection
	465	c======================================================================
[204]	466	print*,'Avant convection'
	467	do it=1,nqmax
	468	WRITE(itn,'(i1)') it
	469	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
	470	enddo
[2]	471
	472	if (convection) then
	473
[204]	474	print*,'Pas de temps dans phytrac : ',pdtphys
[2]	475	DO it=1, nqmax
	476	CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
	477	. pplay, paprs, tr_seri(1,1,it), d_tr_cv)
	478	DO k = 1, nlev
	479	DO i = 1, klon
	480	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k)
	481	ENDDO
	482	ENDDO
[204]	483	c WRITE(itn,'(i1)') it
	484	c CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it='//itn)
[2]	485	ENDDO
[204]	486	c print*,'apres nflxtr'
[2]	487
	488
	489	endif ! convection
[204]	490	c print*,'Apres convection'
	491	c do it=1,nqmax
	492	c WRITE(itn,'(i1)') it
	493	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
	494	c enddo
[2]	495
	496	c======================================================================
	497	c Calcul de l'effet de la couche limite
	498	c======================================================================
[204]	499	c print *,'Avant couchelimite'
	500	c do it=1,nqmax
	501	c WRITE(itn,'(i1)') it
	502	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
	503	c enddo
[2]	504
	505	if (couchelimite) then
	506
	507	DO k = 1, nlev
	508	DO i = 1, klon
	509	delp(i,k) = paprs(i,k)-paprs(i,k+1)
	510	ENDDO
	511	ENDDO
	512
	513	C maf modif pour tenir compte du cas rnpb + traceur
	514	DO it=1, nqmax
[53]	515	c print *,'it',it,clsol(it)
[2]	516	if (clsol(it)) then ! couche limite avec quantite dans le sol calculee
	517	CALL cltracrn(it, pdtphys, yu1, yv1,
	518	e coefh,t_seri,ftsol,pctsrf,
	519	e tr_seri(1,1,it),trs(1,it),
	520	e paprs, pplay, delp,
	521	e masktr(1,it),fshtr(1,it),hsoltr(it),
	522	e tautr(it),vdeptr(it),
	523	e xlat,
	524	s d_tr_cl,d_trs)
	525	DO k = 1, nlev
	526	DO i = 1, klon
	527	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k)
	528	ENDDO
	529	ENDDO
	530	c
	531	c Traceur ds sol
	532	c
	533	DO i = 1, klon
	534	trs(i,it) = trs(i,it) + d_trs(i)
	535	END DO
	536	C
	537	C maf provisoire suppression des prints
	538	C WRITE(itn,'(i1)') it
	539	C CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
	540	else ! couche limite avec flux prescrit
	541	Cmaf provisoire source / traceur a creer
	542	DO i=1, klon
	543	source(i) = 0.0 ! pas de source, pour l'instant
	544	ENDDO
	545	C
	546	CALL cltrac(pdtphys, coefh,t_seri,
	547	s tr_seri(1,1,it), source,
	548	e paprs, pplay, delp,
	549	s d_tr )
	550	DO k = 1, nlev
	551	DO i = 1, klon
	552	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr(i,k)
	553	ENDDO
	554	ENDDO
	555	Cmaf provisoire suppression des prints
	556	Cmaf WRITE(itn,'(i1)') it
	557	cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
	558	endif
	559	ENDDO
	560	c
	561	endif ! couche limite
	562
[204]	563	c print*,'Apres couchelimite'
	564	c do it=1,nqmax
	565	c WRITE(itn,'(i1)') it
	566	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
	567	c enddo
[2]	568
	569	c======================================================================
	570	c Calcul de l'effet du puits radioactif
	571	c======================================================================
	572
	573	C MAF il faudrait faire une modification pour passer dans radiornpb
	574	c si radio=true mais pour l'instant radiornpb propre au cas rnpb
	575	if(rnpb) then
	576	call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
	577	C
	578	DO it=1,nqmax
	579	if(radio(it)) then
	580	DO k = 1, nlev
	581	DO i = 1, klon
	582	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
	583	ENDDO
	584	ENDDO
	585	WRITE(itn,'(i1)') it
	586	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn)
	587	endif
	588	ENDDO
	589	c
	590	endif ! rnpb decroissance radioactive
	591	C
	592	c======================================================================
	593	c Calcul de l'effet de la precipitation
	594	c======================================================================
	595
[204]	596	print*,'LESSIVAGE =',lessivage
[2]	597	IF (lessivage) THEN
	598
[53]	599	c print*,'avant lessivage'
[2]	600
	601	DO it = 1, nqmax
	602	DO k = 1, nlev
	603	DO i = 1, klon
	604	d_tr_lessi_nucl(i,k,it) = 0.0
	605	d_tr_lessi_impa(i,k,it) = 0.0
	606	flestottr(i,k,it) = 0.0
	607	ENDDO
	608	ENDDO
	609	ENDDO
	610	c
	611	c tendance des aerosols nuclees et impactes
	612	c
	613	DO it = 1, nqmax
	614	IF (aerosol(it)) THEN
	615	DO k = 1, nlev
	616	DO i = 1, klon
	617	d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +
	618	s ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
	619	d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +
	620	s ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
	621	ENDDO
	622	ENDDO
	623	ENDIF
	624	ENDDO
	625	c
	626	c Mises a jour des traceurs + calcul des flux de lessivage
	627	c Mise a jour due a l'impaction et a la nucleation
	628	c
[204]	629	c call dump2d(iim,jjm-1,frac_impa(2:klon-1,10),'FRACIMPA')
	630	c call dump2d(iim,jjm-1,frac_nucl(2:klon-1,10),'FRACNUCL')
	631	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3')
[2]	632	DO it = 1, nqmax
[204]	633	c print*,'IT=',it,aerosol(it)
[2]	634	IF (aerosol(it)) THEN
[204]	635	c print*,'IT=',it,' On lessive'
[2]	636	DO k = 1, nlev
	637	DO i = 1, klon
[204]	638	tr_seri(i,k,it)=tr_seri(i,k,it)
	639	s frac_impa(i,k)frac_nucl(i,k)
[2]	640	ENDDO
	641	ENDDO
	642	ENDIF
	643	ENDDO
[204]	644	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3B')
[2]	645	c
	646	c Flux lessivage total
	647	c
	648	DO it = 1, nqmax
	649	DO k = 1, nlev
	650	DO i = 1, klon
	651	flestottr(i,k,it) = flestottr(i,k,it) -
	652	s ( d_tr_lessi_nucl(i,k,it) +
	653	s d_tr_lessi_impa(i,k,it) ) *
	654	s ( paprs(i,k)-paprs(i,k+1) ) /
	655	s (RG * pdtphys)
	656	ENDDO
	657	ENDDO
	658	c
	659	Cmaf suppression provisoire
	660	Cmaf WRITE(itn,'(i1)') it
	661	Cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
	662	ENDDO
	663	c
[53]	664	c print*,'apres lessivage'
[2]	665	ENDIF
	666	Cc
	667	DO k = 1, nlev
	668	DO i = 1, klon
	669	fluxrn(i,k) = flestottr(i,k,1)
	670	fluxpb(i,k) = flestottr(i,k,2)
	671	rn(i,k) = tr_seri(i,k,1)
	672	pb(i,k) = tr_seri(i,k,2)
	673	pbnucl(i,k)=d_tr_lessi_nucl(i,k,2)
	674	pbimpa(i,k)=d_tr_lessi_impa(i,k,2)
	675	ENDDO
	676	ENDDO
	677	itra=itra+1
[204]	678	ndex(1) = 0
[2]	679	DO it=1,nqmax
	680	IF (it.LE.99) THEN
	681	WRITE(str2,'(i2.2)') it
	682	C champs 3D
	683	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,tr_seri(1,1,it),zx_tmp_3d)
	684	CALL histwrite(nid_tra,"tr"//str2,itra,zx_tmp_3d,
[204]	685	. iim(jjm+1)klev,ndex)
	686	c IF (lessivage) THEN
	687	c CALL gr_fi_ecrit(klev,klon,iim,jjm+1,flestottr(1,1,it),zx_tmp_3d)
	688	c CALL histwrite(nid_tra,"fl"//str2,itra,zx_tmp_3d,
	689	c . iim(jjm+1)klev,ndex)
	690	c ENDIF
[2]	691	ELSE
	692	PRINT*, "Trop de traceurs"
	693	CALL abort
	694	ENDIF
	695	ENDDO
	696
[204]	697	goto 888
	698	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pplay,zx_tmp_3d)
	699	CALL histwrite(nid_tra,"pplay",itra,zx_tmp_3d,
	700	. iim(jjm+1)klev,ndex)
	701
	702	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,t_seri,zx_tmp_3d)
	703	CALL histwrite(nid_tra,"t",itra,zx_tmp_3d,
	704	. iim(jjm+1)klev,ndex)
	705	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pmfu,zx_tmp_3d)
	706	CALL histwrite(nid_tra,"mfu",itra,zx_tmp_3d,
	707	. iim(jjm+1)klev,ndex)
	708	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pmfd,zx_tmp_3d)
	709	CALL histwrite(nid_tra,"mfd",itra,zx_tmp_3d,
	710	. iim(jjm+1)klev,ndex)
	711	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pen_u,zx_tmp_3d)
	712	CALL histwrite(nid_tra,"en_u",itra,zx_tmp_3d,
	713	. iim(jjm+1)klev,ndex)
	714	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pen_d,zx_tmp_3d)
	715	CALL histwrite(nid_tra,"en_d",itra,zx_tmp_3d,
	716	. iim(jjm+1)klev,ndex)
	717	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pde_d,zx_tmp_3d)
	718	CALL histwrite(nid_tra,"de_d",itra,zx_tmp_3d,
	719	. iim(jjm+1)klev,ndex)
	720	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,pde_u,zx_tmp_3d)
	721	CALL histwrite(nid_tra,"de_u",itra,zx_tmp_3d,
	722	. iim(jjm+1)klev,ndex)
	723	CALL gr_fi_ecrit(klev,klon,iim,jjm+1,coefh,zx_tmp_3d)
	724	CALL histwrite(nid_tra,"coefh",itra,zx_tmp_3d,
	725	. iim(jjm+1)klev,ndex)
	726
	727	888 continue
	728
	729	c print*,'Sortie phytrac'
	730	c do it=1,nqmax
	731	c WRITE(itn,'(i1)') it
	732	c call diagtracphy(tr_seri(:,:,it),paprs,'Fin Phys '//itn)
	733	c enddo
	734
[199]	735	if (lafin) then
	736	print*, 'c est la fin de la physique'
	737	open (99,file='restarttrac', form='formatted')
	738	do i=1,klon
	739	write(99,*) trs(i,1)
	740	enddo
	741	PRINT*, 'Ecriture du fichier restarttrac'
	742	close(99)
	743	else
	744	print*, 'physique pas fini'
	745	endif
	746
	747
[2]	748	RETURN
	749	END

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

Download in other formats:

Original Format