Context Navigation

phytrac.F @ 1629

Last change on this file since 1629 was 634, checked in by Laurent Fairhead, 20 years ago
Modifications faites à la physique pour la rendre parallele YM Une branche de travail LMDZ4_par_0 a été créée provisoirement afin de tester les modifs pleinement avant leurs inclusions dans le tronc principal LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 28.7 KB

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

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

Download in other formats:

Original Format