Context Navigation

phytrac.F @ 5404

Last change on this file since 5404 was 954, checked in by lsce, 17 years ago

ACo+JG

Ajout du flag aerosol_couple :
false=lecture des sulfates dans un fichier(par defaut) - configuration existante
true=calcul des aerosol par INCA

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

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

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

Download in other formats:

Original Format