Context Navigation

phytrac.F @ 919

Last change on this file since 919 was 868, checked in by Laurent Fairhead, 17 years ago
Preparation du remplacement de la physique utilisee pour l'exercice IPCC_AR4 par la version de la physique avec thermique. On garde le repertoire phylmd pour un petit moment pour que les utilisateurs ne soient pas trop perdus ... phy_IPCC_AR4 = phylmd LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 29.7 KB

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

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

Download in other formats:

Original Format