Context Navigation

phytrac.F @ 448

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

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

Download in other formats:

Original Format