Context Navigation

phytrac.F @ 832

Last change on this file since 832 was 814, checked in by Laurent Fairhead, 17 years ago
Rajout de la physique utilisant les thermiques FH LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 29.0 KB

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

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

Download in other formats:

Original Format