Context Navigation

phytrac.F @ 3840

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

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

Download in other formats:

Original Format