Context Navigation

phytrac.F @ 717

Last change on this file since 717 was 682, checked in by lmdzadmin, 19 years ago
Ajout de bilans pour INCA - AC MAF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 29.0 KB

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

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

Download in other formats:

Original Format