Context Navigation

phytrac.F @ 5440

Last change on this file since 5440 was 759, checked in by lsce, 18 years ago
correction de bug pour le couplage avec le modele inca - ACo
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 28.3 KB

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

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

Download in other formats:

Original Format