Context Navigation

phytrac.F @ 717

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

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

Download in other formats:

Original Format