Context Navigation

phytrac.F @ 3656

Last change on this file since 3656 was 353, checked in by lmdzadmin, 23 years ago

2 changements pour les fichiers histoire:

utilisation de l'entree "rectilineaire" de IOIPSL pour ne plus avoir

lancer ncregular a chaque fois

le calendrier des fichiers histoire est maintenant base sur la date d'initialisation de la simulation plutot que sur la date de depart du

job

en cours

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

File size: 26.1 KB

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

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

Download in other formats:

Original Format