Context Navigation

phytrac.F @ 665

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

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

Download in other formats:

Original Format