Context Navigation

phytrac.F @ 623

Last change on this file since 623 was 620, checked in by lmdzadmin, 20 years ago
Rajout variables KE. MAF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 27.5 KB

Rev	Line
[524]	1	!
	2	! $Header$
	3	!
	4	c
	5	c
[619]	6	SUBROUTINE phytrac (iflag_con,rnpb,nstep,
[524]	7	I julien,gmtime,
	8	I debutphy,lafin,
	9	I nqmax,
	10	I nlon,
	11	I nlev,
	12	I pdtphys,
	13	I u,
	14	I v,
	15	I t_seri,
	16	I paprs,
	17	I pplay,
	18	I pmfu,
	19	I pmfd,
	20	I pen_u,
	21	I pde_u,
	22	I pen_d,
	23	I pde_d,
	24	I coefh,
[541]	25	I fm_therm,entr_therm,
[524]	26	I yu1,
	27	I yv1,
	28	I ftsol,
	29	I pctsrf,
	30	I xlat,
	31	I frac_impa,
	32	I frac_nucl,
	33	I xlon,
	34	I presnivs,
	35	I pphis,
	36	I pphi,
	37	I albsol,
	38	I sh,
	39	I rh,
	40	I cldfra,
	41	I rneb,
	42	I diafra,
	43	I cldliq,
	44	I itop_con,
	45	I ibas_con,
	46	I pmflxr,
	47	I pmflxs,
	48	I prfl,
	49	I psfl,
[620]	50	I da,phi,mp,
	51	I upwd,dnwd,
[616]	52	#ifdef INCA
[524]	53	I flxmass_w,
	54	#endif
	55	O tr_seri)
	56
	57	USE ioipsl
	58
	59	#ifdef INCA
	60	USE sflx
	61	USE chem_tracnm
	62	USE species_names
	63	USE chem_mods
	64	USE pht_tables, ONLY : jrates
	65	USE transport_controls, ONLY : conv_flg, pbl_flg
	66	USE airplane_src, ONLY : ptrop
	67	USE lightning, ONLY : prod_light
	68	#ifdef INCA_AER
	69	USE AEROSOL_MOD, only : ntr,trmx,trnx
[616]	70	USE AEROSOL_DIAG,only : cla,las,tausum,angst,aload,cload,totaerh2o,tau,
	71	$ emiss20,sconc,scavcoef_st,scavcoef_cv
	72	$ ,cload05ss ,cload05bc ,cload05pom ,cload05dust ,cload05so4
	73	$ ,cload125ss ,cload125bc ,cload125pom ,cload125dust ,cload125so4
	74	USE AEROSOL_PROGNOS, ONLY : md,mdw
	75	USE AEROSOL_METEO, only : airm
[524]	76	#endif
[616]	77	#ifdef INCA_NMHC
	78	USE RESISTANCE_DIAGNOSE, ONLY : surf_alb, sol_irrad, surf_temp, surf_wind,
	79	$ aero_resist, lamin_resist, surf_resist
[524]	80	#endif
[616]	81	#endif
[524]	82	IMPLICIT none
	83	c======================================================================
	84	c Auteur(s) FH
	85	c Objet: Moniteur general des tendances traceurs
	86	c
	87	cAA Remarques en vrac:
	88	cAA--------------------
	89	cAA 1/ le call phytrac se fait avec nqmax-2 donc nous avons bien
	90	cAA les vrais traceurs (nbtr) dans phytrac (pas la vapeur ni eau liquide)
	91	cAA 2/ Le choix du radon et du pb se fait juste avec un data
	92	cAA (peu propre). Peut-etre pourrait-on prevoir dans l'avenir
	93	cAA une variable "type de traceur"
	94	c======================================================================
	95	#include "YOMCST.h"
	96	#include "dimensions.h"
	97	#include "dimphy.h"
	98	#include "indicesol.h"
	99	#include "temps.h"
	100	#include "paramet.h"
	101	#include "control.h"
	102	#include "comgeomphy.h"
	103	#include "advtrac.h"
	104	c======================================================================
	105
	106	c Arguments:
	107	c
	108	c EN ENTREE:
	109	c ==========
	110	c
	111	c divers:
	112	c -------
	113	c
	114	integer nlon ! nombre de points horizontaux
	115	integer nlev ! nombre de couches verticales
	116	integer nqmax ! nombre de traceurs auxquels on applique la physique
	117	integer nstep ! appel physique
	118	integer julien !jour julien
	119	integer itop_con(nlon)
	120	integer ibas_con(nlon)
	121	real gmtime
	122	real pdtphys ! pas d'integration pour la physique (seconde)
	123	real t_seri(nlon,nlev) ! temperature
	124	real tr_seri(nlon,nlev,nbtr) ! traceur
	125	real u(klon,klev)
	126	real v(klon,klev)
	127	real sh(nlon,nlev) ! humidite specifique
	128	real rh(nlon,nlev) ! humidite relative
	129	real cldliq(nlon,nlev) ! eau liquide nuageuse
	130	real cldfra(nlon,nlev) ! fraction nuageuse (tous les nuages)
	131	real diafra(nlon,nlev) ! fraction nuageuse (convection ou stratus artificiels)
	132	real rneb(nlon,nlev) ! fraction nuageuse (grande echelle)
	133	real albsol(nlon) ! albedo surface
	134	real paprs(nlon,nlev+1) ! pression pour chaque inter-couche (en Pa)
	135	real ps(nlon) ! pression surface
	136	real pplay(nlon,nlev) ! pression pour le mileu de chaque couche (en Pa)
	137	real pphi(nlon,klev) ! geopotentiel
	138	real pphis(klon)
	139	REAL presnivs(klev)
	140	logical debutphy ! le flag de l'initialisation de la physique
	141	logical lafin ! le flag de la fin de la physique
[541]	142	c Olivia
	143	integer isplit,nsplit
[524]	144	REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1) !--lessivage convection
	145	REAL prfl(klon,klev+1), psfl(klon,klev+1) !--lessivage large-scale
	146
[616]	147	#ifdef INCA
[524]	148	REAL flxmass_w(klon,klev)
	149	#endif
[619]	150	integer iflag_con
[524]	151
	152	cAA Rem : nbtr : nombre de vrais traceurs est defini dans dimphy.h
	153	c
	154	c convection:
	155	c -----------
	156	c
	157	REAL pmfu(nlon,nlev) ! flux de masse dans le panache montant
	158	REAL pmfd(nlon,nlev) ! flux de masse dans le panache descendant
	159	REAL pen_u(nlon,nlev) ! flux entraine dans le panache montant
[541]	160
	161	c
	162	c thermiques:
	163	c -----------
	164	c
	165	real fm_therm(klon,klev+1),entr_therm(klon,klev)
	166	real fm_therm1(klon,klev)
	167	c
[524]	168	REAL pde_u(nlon,nlev) ! flux detraine dans le panache montant
	169	REAL pen_d(nlon,nlev) ! flux entraine dans le panache descendant
	170	REAL pde_d(nlon,nlev) ! flux detraine dans le panache descendant
[619]	171	c KE
	172	real da(nlon,nlev),phi(nlon,nlev,nlev),mp(nlon,nlev)
	173	REAL upwd(nlon,nlev) ! saturated updraft mass flux
	174	REAL dnwd(nlon,nlev) ! saturated downdraft mass flux
	175
[524]	176	c
	177	c Couche limite:
	178	c --------------
	179	c
	180	REAL coefh(nlon,nlev) ! coeff melange CL
	181	REAL yu1(nlon) ! vents au premier niveau
	182	REAL yv1(nlon) ! vents au premier niveau
	183	REAL xlat(nlon) ! latitudes pour chaque point
	184	REAL xlon(nlon) ! longitudes pour chaque point
	185
	186	c
	187	c Lessivage:
	188	c ----------
	189	c
	190	c pour le ON-LINE
	191	c
	192	REAL frac_impa(nlon,nlev) ! fraction d'aerosols impactes
	193	REAL frac_nucl(nlon,nlev) ! fraction d'aerosols nuclees
	194	c
	195	cAA
	196	cAA Arguments necessaires pour les sources et puits de traceur:
	197	cAA ----------------
	198	cAA
	199	real ftsol(nlon,nbsrf) ! Temperature du sol (surf)(Kelvin)
	200	real pctsrf(nlon,nbsrf) ! Pourcentage de sol f(nature du sol)
	201	c abder
	202	real pftsol1(nlon),pftsol2(nlon),pftsol3(nlon),pftsol4(nlon)
	203	real ppsrf1(nlon),ppsrf2(nlon),ppsrf3(nlon),ppsrf4(nlon)
	204	c fin
	205	cAA ----------------------------
	206	cAA VARIABLES LOCALES TRACEURS
	207	cAA ----------------------------
	208	cAA
	209	cAA Sources et puits des traceurs:
	210	cAA ------------------------------
	211	cAA
	212	cAA Pour l'instant seuls les cas du rn et du pb ont ete envisages.
	213
	214	REAL source(klon) ! a voir lorsque le flux est prescrit
	215	cAA
	216	cAA Pour la source de radon et son reservoir de sol
	217	cAA ................................................
	218
	219	REAL trs(klon,nbtr) ! Conc. radon ds le sol
	220	SAVE trs
	221
	222	REAL masktr(klon,nbtr) ! Masque reservoir de sol traceur
	223	c Masque de l'echange avec la surface
	224	c (1 = reservoir) ou (possible => 1 )
	225	SAVE masktr
	226	REAL fshtr(klon,nbtr) ! Flux surfacique dans le reservoir de sol
	227	SAVE fshtr
	228	REAL hsoltr(nbtr) ! Epaisseur equivalente du reservoir de sol
	229	SAVE hsoltr
	230	REAL tautr(nbtr) ! Constante de decroissance radioactive
	231	SAVE tautr
	232	REAL vdeptr(nbtr) ! Vitesse de depot sec dans la couche Brownienne
	233	SAVE vdeptr
	234	REAL scavtr(nbtr) ! Coefficient de lessivage
	235	SAVE scavtr
	236	cAA
	237	CHARACTER*2 itn
	238	C maf ioipsl
	239	CHARACTER*2 str2
	240	INTEGER nhori, nvert
	241	REAL zsto, zout, zjulian
	242	INTEGER nid_tra
	243	SAVE nid_tra
	244	#ifdef INCA_AER
	245	INTEGER nid_tra2,nid_tra3
	246	SAVE nid_tra2,nid_tra3
	247	#endif
	248	c REAL x(klon,klev,nbtr+2) ! traceurs
	249	INTEGER ndex(1)
	250	INTEGER ndex2d(iim(jjm+1)),ndex3d(iim(jjm+1)*klev)
	251	REAL zx_tmp_2d(iim,jjm+1), zx_tmp_3d(iim,jjm+1,klev)
	252	REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1)
	253	c
	254	integer itau_w ! pas de temps ecriture = nstep + itau_phy
	255	c
	256
	257	C
	258	C Variables liees a l'ecriture de la bande histoire : phytrac.nc
	259	c
	260	INTEGER ecrit_tra
	261	SAVE ecrit_tra
	262	logical ok_sync
	263	parameter (ok_sync = .true.)
	264	C
	265	C nature du traceur
	266	c
	267	logical aerosol(nbtr) ! Nature du traceur
	268	c ! aerosol(it) = true => aerosol
	269	c ! aerosol(it) = false => gaz
	270	c ! nat_trac(it) = 1. aerosol
	271	logical clsol(nbtr) ! clsol(it) = true => CL sol calculee
	272	logical radio(nbtr) ! radio(it)=true => decroisssance radioactive
	273	save aerosol,clsol,radio
	274	C
	275	c======================================================================
	276	c
	277	c Declaration des procedures appelees
	278	c
	279	c--modif convection tiedtke
	280	INTEGER i, k, it
	281	INTEGER iq, iiq
	282	REAL delp(klon,klev)
	283	c--end modif
	284	c
	285	c Variables liees a l'ecriture de la bande histoire physique
	286	c
	287	c Variables locales pour effectuer les appels en serie
	288	c----------------------------------------------------
	289	c
	290	REAL d_tr(klon,klev), d_trs(klon) ! tendances de traceurs
[619]	291	REAL d_tr_cl(klon,klev,nbtr) ! tendance de traceurs couche limite
	292	REAL d_tr_cv(klon,klev,nbtr) ! tendance de traceurs conv pour chq traceur
[541]	293	REAL d_tr_th(klon,klev,nbtr) ! la tendance des thermiques
[524]	294	REAL d_tr_dec(klon,klev,nbtr) ! la tendance de la decroissance
	295	c ! radioactive du rn - > pb
	296	REAL d_tr_lessi_impa(klon,klev,nbtr) ! la tendance du lessivage
	297	c ! par impaction
	298	REAL d_tr_lessi_nucl(klon,klev,nbtr) ! la tendance du lessivage
	299	c ! par nucleation
	300	REAL fluxrn(klon,klev)
	301	REAL fluxpb(klon,klev)
	302	REAL pbimpa(klon,klev)
	303	REAL pbnucl(klon,klev)
	304	REAL rn(klon,klev)
	305	REAL pb(klon,klev)
	306	REAL flestottr(klon,klev,nbtr) ! flux de lessivage
	307	c ! dans chaque couche
[541]	308	real zmasse(klon,klev)
	309	real ztra_th(klon,klev)
[524]	310
	311	C
	312	character*20 modname
	313	character*80 abort_message
	314	c
	315	c Controles
	316	c-------------
	317	logical first,couchelimite,convection,lessivage,sorties,
[541]	318	s rnpb,inirnpb,thermiques
	319	save first,couchelimite,convection,lessivage,thermiques,
	320	s sorties,inirnpb
	321	c data first,couchelimite,convection,lessivage,sorties
	322	c s /.true.,.true.,.false.,.true.,.true./
	323	c Olivia
	324	data first,couchelimite,convection,lessivage,
	325	s thermiques,sorties
	326	s /.true.,.true.,.true.,.true.,.true.,.true./
[524]	327
[541]	328
[524]	329	#ifdef INCA
[616]	330	INTEGER :: lastgas
[524]	331	INTEGER :: ncsec
	332
[616]	333	INTEGER :: prt_flag_ts(nbtr)=(/1,1,1
[524]	334	#ifdef INCA_CH4
[616]	335	. ,0,0,1,1,1,1,1,
	336	. 0,1,0,0,0,0,0,1,0,0,
[524]	337	. 0,1,1,1,1,0,1,1,1,0,
	338	. 1,1,1,1,1,1,1,1,1,1,
	339	. 1,0,0
[616]	340	#ifdef INCA_AER
	341	. ,1,1,1,1,0,1,1,1
[524]	342	#endif
[616]	343	#endif
[524]	344	#ifdef INCA_AER
[616]	345	c aerosol tracers
	346	. ,1,0,1,1,1,1,1,1,0,1,
	347	. 0,1,1,1,1,1,0,1,0,1,1,1
[524]	348	#endif
[616]	349	#ifdef INCA_NMHC
	350	. , 1, 1, 1, 1, 1, 1, 1,
	351	. 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
	352	. 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
	353	. 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
	354	. 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
	355	. 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
	356	. 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
	357	. 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
	358	. 1, 1, 1, 1, 1, 1, 1
	359	#endif
[524]	360	. /)
	361
	362
	363	REAL, PARAMETER :: dry_mass = 28.966
	364	REAL, POINTER :: hbuf(:)
	365	REAL, ALLOCATABLE :: obuf(:)
	366	REAL :: calday
	367	REAL :: pdel(klon,klev)
	368	REAL :: dummy(klon,klev) = 0.
	369	#endif
	370	#ifdef INCA_AER
	371	integer la
	372	#endif
	373	c
	374	c======================================================================
	375	modname='phytrac'
	376
	377	ps(:)=paprs(:,1)
	378
	379	if (debutphy) then
	380
	381	ecrit_tra = NINT(86400./pdtphys *ecritphy)
	382	print*,'dans phytrac ',pdtphys,ecritphy,ecrit_tra
	383
	384	if(nbtr.lt.nqmax) then
	385	c print*,'NQMAX=',nqmax
	386	c print*,'NBTR=',nbtr
	387	abort_message='See above'
	388	call abort_gcm(modname,abort_message,1)
	389	endif
	390
	391	inirnpb=rnpb
	392	PRINT*, 'La frequence de sortie traceurs est ', ecrit_tra
	393	C
	394	c=============================================================
	395	c Initialisation des sorties
	396	c=============================================================
	397
	398	#ifdef CPP_IOIPSL
	399	#include "ini_histrac.h"
	400	#endif
	401
	402	c======================================================================
	403	c Initialisation de certaines variables pour le Rn et le Pb
	404	c======================================================================
	405
	406	c Initialisation du traceur dans le sol (couche limite radonique)
	407	c
	408	c print*,'valeur de debut dans phytrac :',debutphy
	409	trs(:,:) = 0.
	410
	411	open (99,file='starttrac',status='old',
	412	. err=999,form='formatted')
	413	read(99,*) (trs(i,1),i=1,klon)
	414	999 close(99)
	415	c print*, 'apres starttrac'
	416
	417	c Initialisation de la fraction d'aerosols lessivee
	418	c
	419	d_tr_lessi_impa(:,:,:) = 0.
	420	d_tr_lessi_nucl(:,:,:) = 0.
	421	c
	422	c Initialisation de la nature des traceurs
	423	c
	424	DO it = 1, nqmax
	425	aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut
	426	radio(it) = .FALSE. ! Par defaut pas de passage par radiornpb
	427	clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit
	428	ENDDO
	429	c
	430	ENDIF ! fin debutphy
	431	c Initialisation du traceur dans le sol (couche limite radonique)
	432	if(inirnpb) THEN
	433	c
	434	radio(1)= .true.
	435	radio(2)= .true.
	436	clsol(1)= .true.
	437	clsol(2)= .true.
	438	aerosol(2) = .TRUE. ! le Pb est un aerosol
	439	c
	440	call initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr
	441	. ,vdeptr,scavtr)
	442	inirnpb=.false.
	443	endif
	444	#ifdef INCA
	445	!======================================================================
	446	! Chimie
	447	!======================================================================
	448
	449	calday = FLOAT(julien) + gmtime
	450	ncsec = NINT (86400.*gmtime)
	451
	452	DO k = 1, nlev
	453	pdel(:,k) = paprs(:,k) - paprs (:,k+1)
	454	END DO
	455
	456	#ifdef INCAINFO
	457	PRINT *, 'CHEMMAIN @ ', calday, ' ... '
	458	DO it = 1, nbtr
	459	PRINT *, solsym(it), MINVAL(tr_seri(:,:,it)),
	460	$ MAXVAL(tr_seri(:,:,it))
	461	END DO
	462	#endif
	463
	464
	465	#ifdef INCA_AER
[616]	466
	467	! Changement Anne 01/04/2005
[524]	468	CALL aerosolmain (tr_seri,
	469	$ pdtphys,
	470	$ pplay,
[616]	471	$ pdel,
[524]	472	$ prfl,
	473	$ pmflxr,
	474	$ psfl,
	475	$ pmflxs,
	476	$ pmfu,
	477	$ itop_con,
	478	$ ibas_con,
	479	$ pphi,
[616]	480	$ airephy, ! paire,
	481	$ nstep,
	482	$ rneb, ! for chimiaq
	483	$ t_seri, ! for chimiaq
	484	$ rh)
	485	! fin changement anne
	486
[524]	487	#endif
	488
	489	CALL chemmain (tr_seri, !mmr
	490	$ nas, !nas
	491	$ nstep, !nstep
	492	$ calday, !calday
	493	$ julien, !ncdate
	494	$ ncsec, !ncsec
	495	$ 1, !lat
	496	$ pdtphys, !delt
	497	$ paprs(1,1), !ps
	498	$ pplay, !pmid
	499	$ pdel, !pdel
[616]	500	$ airephy,
[593]	501	$ pctsrf(1,1),!oro
[524]	502	$ ftsol, !tsurf
	503	$ albsol, !albs
	504	$ pphi, !zma
	505	$ pphis, !phis
	506	$ cldfra, !cldfr
	507	$ rneb, !cldfr_st
	508	$ diafra, !cldfr_cv
	509	$ itop_con, !cldtop
	510	$ ibas_con, !cldbot
	511	$ cldliq, !cwat
	512	$ prfl, !flxrst
	513	$ pmflxr, !flxrcv
	514	$ psfl, !flxsst
	515	$ pmflxs, !flxscv
	516	$ pmfu, !flxupd
	517	$ flxmass_w, !flxmass_w
	518	$ t_seri, !tfld
	519	$ sh, !sh
	520	$ rh, !rh
	521	$ .false., !wrhstts
	522	$ hbuf, !hbuf
	523	$ obuf, !obuf
	524	$ iip1, !nx
	525	$ jjp1) !ny
	526	#ifdef INCAINFO
	527	PRINT *, 'OK.'
	528	DO it = 1, nbtr
	529	PRINT *, solsym(it), MINVAL(tr_seri(:,:,it)),
	530	$ MAXVAL(tr_seri(:,:,it))
	531	END DO
	532	#endif
	533	#else
	534
	535	c Abder
[619]	536	ctestmaf if(nqmax.gt.2) aerosol(3)=.true.
[524]	537
	538	do i=1,nlon
	539	pftsol1(i) = ftsol(i,1)
	540	pftsol2(i) = ftsol(i,2)
	541	pftsol3(i) = ftsol(i,3)
	542	pftsol4(i) = ftsol(i,4)
	543
	544	ppsrf1(i) = pctsrf(i,1)
	545	ppsrf2(i) = pctsrf(i,2)
	546	ppsrf3(i) = pctsrf(i,3)
	547	ppsrf4(i) = pctsrf(i,4)
	548
	549	enddo
	550	c Abder
	551	#endif
	552	c======================================================================
	553	c Calcul de l'effet de la convection
	554	c======================================================================
	555	c print*,'Avant convection'
[619]	556	c do it=1,nqmax
	557	c WRITE(itn,'(i2)') it
[524]	558	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
[619]	559	c enddo
[524]	560
	561	if (convection) then
	562
	563	c print*,'Pas de temps dans phytrac : ',pdtphys
	564	DO it=1, nqmax
	565	#ifdef INCA
	566	IF ( conv_flg(it) == 0 ) CYCLE
	567	#endif
[619]	568	if (iflag_con.eq.2) then
	569	c tiedke
[524]	570	CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
[619]	571	. pplay, paprs, tr_seri(1,1,it), d_tr_cv(1,1,it))
	572	else if (iflag_con.eq.3) then
	573	c KE
	574	call cvltr(pdtphys, da, phi, mp, paprs,pplay, tr_seri(1,1,it),
	575	. upwd,dnwd,d_tr_cv(1,1,it))
	576	endif
	577
[524]	578	DO k = 1, nlev
	579	DO i = 1, klon
[619]	580	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it)
[524]	581	ENDDO
	582	ENDDO
	583	#ifdef INCA
	584	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '
	585	. //solsym(it))
	586	#else
	587	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '//itn)
	588	#endif
	589	ENDDO
	590	c print*,'apres nflxtr'
	591
	592	endif ! convection
	593	c print*,'Apres convection'
	594	c do it=1,nqmax
	595	c WRITE(itn,'(i1)') it
	596	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
	597	c enddo
	598
[541]	599
[524]	600	c======================================================================
[541]	601	c Calcul de l'effet des thermiques
	602	c======================================================================
	603
	604	do k=1,klev
	605	do i=1,klon
	606	zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
	607	enddo
	608	enddo
	609
	610	c print*,'masse dans ph ',zmasse
	611	do it=1,nqmax
	612	do k=1,klev
	613	do i=1,klon
	614	d_tr_th(i,k,it)=0.
	615	tr_seri(i,k,it)=max(tr_seri(i,k,it),0.)
	616	tr_seri(i,k,it)=min(tr_seri(i,k,it),1.e10)
	617	enddo
	618	enddo
	619	enddo
	620
	621	if (thermiques) then
	622	print*,'calcul de leffet des thermiques'
	623	nsplit=10
	624	DO it=1, nqmax
	625	c WRITE(itn,'(i1)') it
	626	c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'conv it='//itn)
	627	c print*,'avant dqthermiquesretro'
	628	c call dump2d(iim,jjm-1,tr_seri(2,1,1),'TR_SERI ')
	629
	630	do isplit=1,nsplit
	631	c Abderr 25 11 02
	632	C Thermiques
	633	c print*,'Avant dans phytrac',avant
	634	call dqthermcell(klon,klev,pdtphys/nsplit
	635	. ,fm_therm,entr_therm,zmasse
	636	. ,tr_seri(1:klon,1:klev,it),d_tr,ztra_th)
	637
	638	do k=1,klev
	639	do i=1,klon
	640	d_tr(i,k)=pdtphys*d_tr(i,k)/nsplit
	641	d_tr_th(i,k,it)=d_tr_th(i,k,it)+d_tr(i,k)
	642	tr_seri(i,k,it)=max(tr_seri(i,k,it)+d_tr(i,k),0.)
	643	enddo
	644	enddo
	645	enddo ! nsplit
	646	print*,'apres thermiques'
	647	c call dump2d(iim,jjm-1,d_tr_th(1,1,1),'d_tr_th ')
	648	c do k=1,klev
	649	c print*,'d_tr_th(',k,')=',tr_seri(280,k,1)
	650	c enddo
	651
	652	c WRITE(itn,'(i1)') it
	653	c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'therm it='//itn)
	654	ENDDO ! it
	655	endif ! Thermiques
	656	c print*,'ATTENTION: sdans thermniques'
	657
	658	c======================================================================
[524]	659	c Calcul de l'effet de la couche limite
	660	c======================================================================
	661	c print *,'Avant couchelimite'
	662	c do it=1,nqmax
	663	c WRITE(itn,'(i1)') it
	664	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
	665	c enddo
	666
	667	if (couchelimite) then
	668
	669	DO k = 1, nlev
	670	DO i = 1, klon
	671	delp(i,k) = paprs(i,k)-paprs(i,k+1)
	672	ENDDO
	673	ENDDO
	674
	675	C maf modif pour tenir compte du cas rnpb + traceur
	676	DO it=1, nqmax
	677	#ifdef INCA
	678	IF ( pbl_flg(it) == 0 ) CYCLE
	679	#endif
	680	c print *,'it',it,clsol(it)
	681	if (clsol(it)) then ! couche limite avec quantite dans le sol calculee
	682	CALL cltracrn(it, pdtphys, yu1, yv1,
	683	e coefh,t_seri,ftsol,pctsrf,
	684	e tr_seri(1,1,it),trs(1,it),
	685	e paprs, pplay, delp,
	686	e masktr(1,it),fshtr(1,it),hsoltr(it),
	687	e tautr(it),vdeptr(it),
	688	e xlat,
[619]	689	s d_tr_cl(1,1,it),d_trs)
[524]	690	DO k = 1, nlev
	691	DO i = 1, klon
[619]	692	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
[524]	693	ENDDO
	694	ENDDO
	695	c
	696	c Traceur ds sol
	697	c
	698	DO i = 1, klon
	699	trs(i,it) = trs(i,it) + d_trs(i)
	700	END DO
	701	C
	702	C maf provisoire suppression des prints
	703	C WRITE(itn,'(i1)') it
	704	C CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
	705	else ! couche limite avec flux prescrit
	706	#ifdef INCA
	707	DO k = 1, klon
	708	source(k) = eflux(k,it)-dflux(k,it)
	709	END DO
	710	#else
	711
	712	Cmaf provisoire source / traceur a creer
	713	DO i=1, klon
	714	source(i) = 0.0 ! pas de source, pour l'instant
	715	ENDDO
	716	C
	717	#endif
	718	CALL cltrac(pdtphys, coefh,t_seri,
	719	s tr_seri(1,1,it), source,
	720	e paprs, pplay, delp,
[619]	721	s d_tr_cl(1,1,it))
[524]	722	DO k = 1, nlev
	723	DO i = 1, klon
[619]	724	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
[524]	725	ENDDO
	726	ENDDO
	727	Cmaf WRITE(itn,'(i1)') it
	728	cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
	729	endif
	730	ENDDO
	731	c
	732	endif ! couche limite
	733
	734	c print*,'Apres couchelimite'
	735	c do it=1,nqmax
	736	c WRITE(itn,'(i1)') it
	737	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
	738	c enddo
	739
	740	c======================================================================
	741	c Calcul de l'effet du puits radioactif
	742	c======================================================================
	743
	744	C MAF il faudrait faire une modification pour passer dans radiornpb
	745	c si radio=true mais pour l'instant radiornpb propre au cas rnpb
	746	if(rnpb) then
	747	print *, 'decroissance radiactive activee'
	748	call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
	749	C
	750	DO it=1,nqmax
	751	if(radio(it)) then
	752	DO k = 1, nlev
	753	DO i = 1, klon
	754	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
	755	ENDDO
	756	ENDDO
	757	WRITE(itn,'(i1)') it
	758	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn)
	759	endif
	760	ENDDO
	761	c
	762	endif ! rnpb decroissance radioactive
	763	C
	764	c======================================================================
	765	c Calcul de l'effet de la precipitation
	766	c======================================================================
	767
	768	c print*,'LESSIVAGE =',lessivage
	769	IF (lessivage) THEN
	770
	771	c print*,'avant lessivage'
	772
	773	d_tr_lessi_nucl(:,:,:) = 0.
	774	d_tr_lessi_impa(:,:,:) = 0.
	775	flestottr(:,:,:) = 0.
	776	c
	777	c tendance des aerosols nuclees et impactes
	778	c
	779	DO it = 1, nqmax
	780	IF (aerosol(it)) THEN
	781	DO k = 1, nlev
	782	DO i = 1, klon
	783	d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +
	784	s ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
	785	d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +
	786	s ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
	787	ENDDO
	788	ENDDO
	789	ENDIF
	790	ENDDO
	791	c
	792	c Mises a jour des traceurs + calcul des flux de lessivage
	793	c Mise a jour due a l'impaction et a la nucleation
	794	c
	795	c call dump2d(iim,jjm-1,frac_impa(2:klon-1,10),'FRACIMPA')
	796	c call dump2d(iim,jjm-1,frac_nucl(2:klon-1,10),'FRACNUCL')
	797	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3')
	798	DO it = 1, nqmax
	799	c print*,'IT=',it,aerosol(it)
	800	IF (aerosol(it)) THEN
	801	c print*,'IT=',it,' On lessive'
	802	DO k = 1, nlev
	803	DO i = 1, klon
	804	tr_seri(i,k,it)=tr_seri(i,k,it)
	805	s frac_impa(i,k)frac_nucl(i,k)
	806	ENDDO
	807	ENDDO
	808	ENDIF
	809	ENDDO
	810	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3B')
	811	c
	812	c Flux lessivage total
	813	c
	814	DO it = 1, nqmax
	815	DO k = 1, nlev
	816	DO i = 1, klon
	817	flestottr(i,k,it) = flestottr(i,k,it) -
	818	s ( d_tr_lessi_nucl(i,k,it) +
	819	s d_tr_lessi_impa(i,k,it) ) *
	820	s ( paprs(i,k)-paprs(i,k+1) ) /
	821	s (RG * pdtphys)
	822	ENDDO
	823	ENDDO
	824	c
	825	Cmaf WRITE(itn,'(i1)') it
	826	Cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
	827	ENDDO
	828	c
	829	c print*,'apres lessivage'
	830	ENDIF
	831	Cc
	832	DO k = 1, nlev
	833	DO i = 1, klon
	834	fluxrn(i,k) = flestottr(i,k,1)
	835	fluxpb(i,k) = flestottr(i,k,2)
	836	rn(i,k) = tr_seri(i,k,1)
	837	pb(i,k) = tr_seri(i,k,2)
	838	pbnucl(i,k)=d_tr_lessi_nucl(i,k,2)
	839	pbimpa(i,k)=d_tr_lessi_impa(i,k,2)
	840	ENDDO
	841	ENDDO
	842
	843	c=============================================================
	844	c Ecriture des sorties
	845	c=============================================================
	846
	847	#ifdef CPP_IOIPSL
	848	#include "write_histrac.h"
	849	#endif
	850
	851	c=============================================================
	852
	853	if (lafin) then
	854	print*, 'c est la fin de la physique'
	855	open (99,file='restarttrac', form='formatted')
	856	do i=1,klon
	857	write(99,*) trs(i,1)
	858	enddo
	859	PRINT*, 'Ecriture du fichier restarttrac'
	860	close(99)
	861	else
	862	print*, 'physique pas fini'
	863	endif
	864
	865
	866	RETURN
	867	END

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

Download in other formats:

Original Format