Context Navigation

phytrac.F @ 1167

Last change on this file since 1167 was 1150, checked in by jghattas, 17 years ago

Ajoute possiblilite de forcer avec d'autre aersols avec le meme principe que pour les so4.

Par default la convergence est rempu avec flag_aersosol=1 + ok_ade=ok_aie=y. Exactement les memes resultats peuvent etre retrouver avec new_aod=n. Les resultats sont exactement les memes sans prise en compte aucun aerosol, avec ok_ade=ok_aie=n.

flag_aerosol indique quel aersosol ou combinasion a utiliser (=1 uniquement les SO4 comme avant)
les fichiers d'entree so4.runXXXX.cdf change du nom pour majescule SO4.runXXXX.cdf

readsulfate change du nom pour readaerosol qui trait tous les aerosols
radlwsw change du nom pour radlwsw_aero.
aeropt_2bands.F90 et aeropt_5wv.F90 correspond a un reecriture de aeropt.F (premier et deuxieme moitie respectivement)
aeropt.F est gardé pour retrouver la convergence si demande avec new_aod=false (=true default)
sw_aero est un version evolue de sw_lmdar4 (dans fichier radiation_ar4.F)

ACo + JG

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

File size: 26.4 KB

Rev	Line
[524]	1	!
	2	c
	3	c
[658]	4	SUBROUTINE phytrac (rnpb,
[625]	5	I nstep,
	6	I julien,
	7	I gmtime,
	8	I debutphy,
	9	I lafin,
[524]	10	I nlon,
	11	I nlev,
	12	I pdtphys,
	13	I u,
	14	I v,
	15	I t_seri,
	16	I paprs,
	17	I pplay,
	18	I pmfu,
	19	I pmfd,
	20	I pen_u,
	21	I pde_u,
	22	I pen_d,
	23	I pde_d,
[1067]	24	I cdragh,
[524]	25	I coefh,
[625]	26	I fm_therm,
	27	I entr_therm,
[524]	28	I yu1,
	29	I yv1,
	30	I ftsol,
	31	I pctsrf,
	32	I xlat,
	33	I frac_impa,
	34	I frac_nucl,
	35	I xlon,
	36	I presnivs,
	37	I pphis,
	38	I pphi,
	39	I albsol,
	40	I sh,
	41	I rh,
	42	I cldfra,
	43	I rneb,
	44	I diafra,
	45	I cldliq,
	46	I itop_con,
	47	I ibas_con,
	48	I pmflxr,
	49	I pmflxs,
	50	I prfl,
	51	I psfl,
[625]	52	I da,
	53	I phi,
	54	I mp,
	55	I upwd,
	56	I dnwd,
[955]	57	I aerosol_couple,
[524]	58	I flxmass_w,
[1150]	59	I tau_aero,
	60	I piz_aero,
	61	I cg_aero,
[959]	62	I ccm,
	63	I rfname,
[524]	64	O tr_seri)
	65
	66	USE ioipsl
[766]	67	USE dimphy
[1114]	68	USE infotrac
[776]	69	USE mod_grid_phy_lmdz
	70	USE mod_phys_lmdz_para
[766]	71	USE comgeomphy
	72	USE iophy
	73	USE vampir
[524]	74
	75	IMPLICIT none
	76	c======================================================================
	77	c Auteur(s) FH
	78	c Objet: Moniteur general des tendances traceurs
	79	c
	80	cAA Remarques en vrac:
	81	cAA--------------------
	82	cAA 2/ Le choix du radon et du pb se fait juste avec un data
	83	cAA (peu propre). Peut-etre pourrait-on prevoir dans l'avenir
	84	cAA une variable "type de traceur"
	85	c======================================================================
	86	#include "YOMCST.h"
	87	#include "dimensions.h"
	88	#include "indicesol.h"
[658]	89	#include "clesphys.h"
[524]	90	#include "temps.h"
	91	#include "paramet.h"
	92	#include "control.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 nstep ! appel physique
	107	integer julien !jour julien
	108	integer itop_con(nlon)
	109	integer ibas_con(nlon)
	110	real gmtime
	111	real pdtphys ! pas d'integration pour la physique (seconde)
	112	real t_seri(nlon,nlev) ! temperature
	113	real tr_seri(nlon,nlev,nbtr) ! traceur
	114	real u(klon,klev)
	115	real v(klon,klev)
	116	real sh(nlon,nlev) ! humidite specifique
	117	real rh(nlon,nlev) ! humidite relative
	118	real cldliq(nlon,nlev) ! eau liquide nuageuse
	119	real cldfra(nlon,nlev) ! fraction nuageuse (tous les nuages)
	120	real diafra(nlon,nlev) ! fraction nuageuse (convection ou stratus artificiels)
	121	real rneb(nlon,nlev) ! fraction nuageuse (grande echelle)
	122	real albsol(nlon) ! albedo surface
	123	real paprs(nlon,nlev+1) ! pression pour chaque inter-couche (en Pa)
	124	real ps(nlon) ! pression surface
	125	real pplay(nlon,nlev) ! pression pour le mileu de chaque couche (en Pa)
	126	real pphi(nlon,klev) ! geopotentiel
	127	real pphis(klon)
	128	REAL presnivs(klev)
	129	logical debutphy ! le flag de l'initialisation de la physique
	130	logical lafin ! le flag de la fin de la physique
[541]	131	c Olivia
[677]	132	integer nsplit
[524]	133	REAL pmflxr(klon,klev+1), pmflxs(klon,klev+1) !--lessivage convection
	134	REAL prfl(klon,klev+1), psfl(klon,klev+1) !--lessivage large-scale
[955]	135	LOGICAL aerosol_couple
[959]	136
[524]	137	REAL flxmass_w(klon,klev)
[1114]	138	CHARACTER(len=8) :: solsym(nbtr)
[864]	139	integer la
[1150]	140	REAL :: tau_aero(klon,klev,9,2)
	141	REAL :: piz_aero(klon,klev,9,2)
	142	REAL :: cg_aero(klon,klev,9,2)
[864]	143	character*4 :: rfname(9)
	144	REAL :: ccm(klon,klev,2)
[524]	145	c
	146	c convection:
	147	c -----------
	148	c
	149	REAL pmfu(nlon,nlev) ! flux de masse dans le panache montant
	150	REAL pmfd(nlon,nlev) ! flux de masse dans le panache descendant
	151	REAL pen_u(nlon,nlev) ! flux entraine dans le panache montant
[541]	152
	153	c
	154	c thermiques:
	155	c -----------
	156	c
	157	real fm_therm(klon,klev+1),entr_therm(klon,klev)
	158	real fm_therm1(klon,klev)
	159	c
[524]	160	REAL pde_u(nlon,nlev) ! flux detraine dans le panache montant
	161	REAL pen_d(nlon,nlev) ! flux entraine dans le panache descendant
	162	REAL pde_d(nlon,nlev) ! flux detraine dans le panache descendant
[619]	163	c KE
	164	real da(nlon,nlev),phi(nlon,nlev,nlev),mp(nlon,nlev)
	165	REAL upwd(nlon,nlev) ! saturated updraft mass flux
	166	REAL dnwd(nlon,nlev) ! saturated downdraft mass flux
	167
[524]	168	c
	169	c Couche limite:
	170	c --------------
	171	c
[1067]	172	REAL cdragh(nlon,nlev)! coeff drag pour T et Q
[524]	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
[1114]	207	REAL source(klon,nbtr) ! a voir lorsque le flux est prescrit
[524]	208	cAA
	209	cAA Pour la source de radon et son reservoir de sol
	210	cAA ................................................
	211
[766]	212	REAL,save,allocatable :: trs(:,:) ! Conc. radon ds le sol
	213	c$OMP THREADPRIVATE(trs)
[776]	214	REAL :: trs_tmp(klon_glo)
[766]	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 )
[766]	219	c$OMP THREADPRIVATE(masktr)
	220	REAL,save,allocatable :: fshtr(:,:) ! Flux surfacique dans le reservoir de sol
	221	c$OMP THREADPRIVATE(fshtr)
	222	REAL,save,allocatable :: hsoltr(:) ! Epaisseur equivalente du reservoir de sol
	223	c$OMP THREADPRIVATE(hsoltr)
	224	REAL,save,allocatable :: tautr(:) ! Constante de decroissance radioactive
	225	c$OMP THREADPRIVATE(tautr)
	226	REAL,save,allocatable :: vdeptr(:) ! Vitesse de depot sec dans la couche Brownienne
	227	c$OMP THREADPRIVATE(vdeptr)
	228	REAL,save,allocatable :: scavtr(:) ! Coefficient de lessivage
	229	c$OMP THREADPRIVATE(scavtr)
[524]	230	cAA
	231	CHARACTER*2 itn
	232	C maf ioipsl
	233	CHARACTER*2 str2
	234	INTEGER nhori, nvert
	235	REAL zsto, zout, zjulian
[959]	236	INTEGER, SAVE :: nid_tra
[766]	237	c$OMP THREADPRIVATE(nid_tra)
[959]	238	INTEGER, SAVE :: nid_tra2,nid_tra3
[766]	239	c$OMP THREADPRIVATE(nid_tra2,nid_tra3)
[524]	240	INTEGER ndex(1)
	241	INTEGER ndex2d(iim(jjm+1)),ndex3d(iim(jjm+1)*klev)
[766]	242	REAL zx_tmp_fi2d(klon) ! variable temporaire grille physique
	243	REAL zx_tmp_fi3d(klon,klev) ! variable temporaire pour champs 3D
[524]	244	REAL zx_tmp_2d(iim,jjm+1), zx_tmp_3d(iim,jjm+1,klev)
	245	REAL zx_lon(iim,jjm+1), zx_lat(iim,jjm+1)
	246	c
	247	integer itau_w ! pas de temps ecriture = nstep + itau_phy
	248	c
	249	logical ok_sync
	250	parameter (ok_sync = .true.)
	251	C
	252	C nature du traceur
	253	c
[766]	254	logical,save,allocatable :: aerosol(:) ! Nature du traceur
[524]	255	c ! aerosol(it) = true => aerosol
	256	c ! aerosol(it) = false => gaz
	257	c ! nat_trac(it) = 1. aerosol
[766]	258	logical,save,allocatable :: clsol(:) ! clsol(it) = true => CL sol calculee
	259	logical,save,allocatable :: radio(:) ! radio(it)=true => decroisssance radioactive
	260	c$OMP THREADPRIVATE(aerosol,clsol,radio)
[524]	261	C
	262	c======================================================================
	263	c
	264	c Declaration des procedures appelees
	265	c
	266	c--modif convection tiedtke
	267	INTEGER i, k, it
	268	INTEGER iq, iiq
	269	REAL delp(klon,klev)
	270	c--end modif
	271	c
	272	c Variables liees a l'ecriture de la bande histoire physique
	273	c
	274	c Variables locales pour effectuer les appels en serie
	275	c----------------------------------------------------
	276	c
	277	REAL d_tr(klon,klev), d_trs(klon) ! tendances de traceurs
[619]	278	REAL d_tr_cl(klon,klev,nbtr) ! tendance de traceurs couche limite
	279	REAL d_tr_cv(klon,klev,nbtr) ! tendance de traceurs conv pour chq traceur
[541]	280	REAL d_tr_th(klon,klev,nbtr) ! la tendance des thermiques
[524]	281	REAL d_tr_dec(klon,klev,nbtr) ! la tendance de la decroissance
	282	c ! radioactive du rn - > pb
	283	REAL d_tr_lessi_impa(klon,klev,nbtr) ! la tendance du lessivage
	284	c ! par impaction
	285	REAL d_tr_lessi_nucl(klon,klev,nbtr) ! la tendance du lessivage
	286	c ! par nucleation
	287	REAL fluxrn(klon,klev)
	288	REAL fluxpb(klon,klev)
	289	REAL pbimpa(klon,klev)
	290	REAL pbnucl(klon,klev)
	291	REAL rn(klon,klev)
	292	REAL pb(klon,klev)
	293	REAL flestottr(klon,klev,nbtr) ! flux de lessivage
	294	c ! dans chaque couche
[541]	295	real zmasse(klon,klev)
	296	real ztra_th(klon,klev)
[524]	297
	298	C
	299	character*20 modname
	300	character*80 abort_message
	301	c
	302	c Controles
	303	c-------------
	304	logical first,couchelimite,convection,lessivage,sorties,
[662]	305	s rnpb,inirnpb
	306	save first,couchelimite,convection,lessivage,
[541]	307	s sorties,inirnpb
[766]	308	c$OMP THREADPRIVATE(first,couchelimite,convection,lessivage,
	309	c$OMP+ sorties,inirnpb)
[541]	310	c data first,couchelimite,convection,lessivage,sorties
	311	c s /.true.,.true.,.false.,.true.,.true./
	312	c Olivia
	313	data first,couchelimite,convection,lessivage,
[662]	314	s sorties
	315	s /.true.,.true.,.true.,.true.,.true./
[524]	316
[959]	317	! Variables needed for configuration with INCA
[616]	318	INTEGER :: lastgas
[524]	319	INTEGER :: ncsec
	320
[766]	321	REAL, PARAMETER :: dry_mass = 28.966
	322	REAL, POINTER :: hbuf(:)
	323	REAL, ALLOCATABLE :: obuf(:)
	324	REAL :: calday
	325	REAL :: pdel(klon,klev)
	326	c
	327	c======================================================================
	328
[524]	329	modname='phytrac'
	330
	331	ps(:)=paprs(:,1)
	332
	333	if (debutphy) then
[766]	334	allocate( trs(klon,nbtr) )
	335	allocate( masktr(klon,nbtr))
	336	allocate( fshtr(klon,nbtr) )
	337	allocate( hsoltr(nbtr))
	338	allocate( tautr(nbtr))
	339	allocate( vdeptr(nbtr))
	340	allocate( scavtr(nbtr))
	341	allocate( aerosol(nbtr))
	342	allocate( clsol(nbtr))
	343	allocate( radio(nbtr))
[524]	344
[766]	345
[956]	346	! FH 2008/05/09 correction de la frequence d'ecriture des traceurs
	347	! ecrit_tra = FLOAT(NINT(86400./pdtphys *ecritphy))
	348	print*,'dans phytrac ',pdtphys,ecrit_tra
[524]	349
	350	inirnpb=rnpb
	351	PRINT*, 'La frequence de sortie traceurs est ', ecrit_tra
	352	C
	353	c=============================================================
	354	c Initialisation des sorties
	355	c=============================================================
	356
	357	#ifdef CPP_IOIPSL
	358	#include "ini_histrac.h"
	359	#endif
	360
	361	c======================================================================
	362	c Initialisation de certaines variables pour le Rn et le Pb
	363	c======================================================================
	364
	365	c Initialisation du traceur dans le sol (couche limite radonique)
	366	c
	367	c print*,'valeur de debut dans phytrac :',debutphy
	368	trs(:,:) = 0.
[766]	369	c$OMP MASTER
[776]	370	if (is_mpi_root) then
[766]	371	trs_tmp(:)=0.
	372	open (99,file='starttrac',status='old',
[524]	373	. err=999,form='formatted')
[776]	374	read(99,*) (trs_tmp(i),i=1,klon_glo)
[524]	375	999 close(99)
[766]	376	endif
	377	c$OMP END MASTER
[776]	378	call Scatter(trs_tmp,trs(:,1))
	379
[524]	380	c print*, 'apres starttrac'
	381
	382	c Initialisation de la fraction d'aerosols lessivee
	383	c
	384	d_tr_lessi_impa(:,:,:) = 0.
	385	d_tr_lessi_nucl(:,:,:) = 0.
	386	c
	387	c Initialisation de la nature des traceurs
	388	c
[1114]	389	DO it = 1, nbtr
[524]	390	aerosol(it) = .FALSE. ! Tous les traceurs sont des gaz par defaut
	391	radio(it) = .FALSE. ! Par defaut pas de passage par radiornpb
	392	clsol(it) = .FALSE. ! Par defaut couche limite avec flux prescrit
	393	ENDDO
	394	c
	395	ENDIF ! fin debutphy
	396	c Initialisation du traceur dans le sol (couche limite radonique)
	397	if(inirnpb) THEN
	398	c
	399	radio(1)= .true.
	400	radio(2)= .true.
	401	clsol(1)= .true.
	402	clsol(2)= .true.
	403	aerosol(2) = .TRUE. ! le Pb est un aerosol
	404	c
	405	call initrrnpb (ftsol,pctsrf,masktr,fshtr,hsoltr,tautr
	406	. ,vdeptr,scavtr)
	407	inirnpb=.false.
	408	endif
[959]	409
	410
	411	IF (config_inca == 'none') THEN
	412	DO i=1,nlon
	413	pftsol1(i) = ftsol(i,1)
	414	pftsol2(i) = ftsol(i,2)
	415	pftsol3(i) = ftsol(i,3)
	416	pftsol4(i) = ftsol(i,4)
	417
	418	ppsrf1(i) = pctsrf(i,1)
	419	ppsrf2(i) = pctsrf(i,2)
	420	ppsrf3(i) = pctsrf(i,3)
	421	ppsrf4(i) = pctsrf(i,4)
	422
	423	ENDDO
	424
	425	ELSE ! config_inca /=none
[524]	426	#ifdef INCA
[766]	427	call VTe(VTphysiq)
	428	call VTb(VTinca)
[959]	429
[524]	430	!======================================================================
	431	! Chimie
	432	!======================================================================
	433
	434	calday = FLOAT(julien) + gmtime
	435	ncsec = NINT (86400.*gmtime)
	436
	437	DO k = 1, nlev
	438	pdel(:,k) = paprs(:,k) - paprs (:,k+1)
	439	END DO
	440
[959]	441	IF (config_inca == 'aero') THEN
	442	CALL aerosolmain (aerosol_couple,
[955]	443	$ tr_seri,
[524]	444	$ pdtphys,
	445	$ pplay,
[616]	446	$ pdel,
[524]	447	$ prfl,
	448	$ pmflxr,
	449	$ psfl,
	450	$ pmflxs,
	451	$ pmfu,
	452	$ itop_con,
	453	$ ibas_con,
	454	$ pphi,
[616]	455	$ airephy, ! paire,
	456	$ nstep,
	457	$ rneb, ! for chimiaq
	458	$ t_seri, ! for chimiaq
[766]	459	$ rh,
[1150]	460	$ tau_aero,
	461	$ piz_aero,
	462	$ cg_aero,
[864]	463	$ rfname,
	464	$ ccm,
[766]	465	$ lafin)
[959]	466	END IF
[524]	467
	468	CALL chemmain (tr_seri, !mmr
	469	$ nstep, !nstep
	470	$ calday, !calday
	471	$ julien, !ncdate
	472	$ ncsec, !ncsec
	473	$ 1, !lat
	474	$ pdtphys, !delt
	475	$ paprs(1,1), !ps
	476	$ pplay, !pmid
	477	$ pdel, !pdel
[616]	478	$ airephy,
[593]	479	$ pctsrf(1,1),!oro
[524]	480	$ ftsol, !tsurf
	481	$ albsol, !albs
	482	$ pphi, !zma
	483	$ pphis, !phis
	484	$ cldfra, !cldfr
	485	$ rneb, !cldfr_st
	486	$ diafra, !cldfr_cv
	487	$ itop_con, !cldtop
	488	$ ibas_con, !cldbot
	489	$ cldliq, !cwat
	490	$ prfl, !flxrst
	491	$ pmflxr, !flxrcv
	492	$ psfl, !flxsst
	493	$ pmflxs, !flxscv
	494	$ pmfu, !flxupd
	495	$ flxmass_w, !flxmass_w
	496	$ t_seri, !tfld
	497	$ sh, !sh
	498	$ rh, !rh
	499	$ iip1, !nx
[766]	500	$ jjp1, !ny
	501	$ source,
	502	$ solsym)
[682]	503
	504
[766]	505	call VTe(VTinca)
	506	call VTb(VTphysiq)
[524]	507	#endif
[959]	508	END IF ! config_inca
[524]	509	c======================================================================
	510	c Calcul de l'effet de la convection
	511	c======================================================================
	512	c print*,'Avant convection'
[1114]	513	do it=1,nbtr
[928]	514	WRITE(itn,'(i2)') it
[524]	515	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
[928]	516	enddo
[524]	517
	518	if (convection) then
	519
	520	c print*,'Pas de temps dans phytrac : ',pdtphys
[1114]	521	DO it=1, nbtr
[959]	522
[1114]	523	IF ( config_inca/='none') THEN
	524	IF ( conv_flg(it) == 0 ) CYCLE
	525	END IF
[959]	526
[928]	527	if (iflag_con.lt.2) then
	528	d_tr_cv=0.
	529	else if (iflag_con.eq.2) then
[619]	530	c tiedke
[524]	531	CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
[619]	532	. pplay, paprs, tr_seri(1,1,it), d_tr_cv(1,1,it))
[928]	533	else
[619]	534	c KE
	535	call cvltr(pdtphys, da, phi, mp, paprs,pplay, tr_seri(1,1,it),
	536	. upwd,dnwd,d_tr_cv(1,1,it))
	537	endif
	538
[928]	539	DO k = 1, nlev
	540	DO i = 1, klon
[619]	541	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it)
[928]	542	ENDDO
	543	ENDDO
[959]	544
	545	IF (config_inca == 'none') THEN
	546	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '//itn)
	547	ELSE
	548	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '
	549	. //solsym(it))
	550	END IF
	551
[524]	552	ENDDO
	553
	554	endif ! convection
	555	c print*,'Apres convection'
[1114]	556	c do it=1,nbtr
[524]	557	c WRITE(itn,'(i1)') it
	558	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
	559	c enddo
	560
[541]	561
[524]	562	c======================================================================
[541]	563	c Calcul de l'effet des thermiques
	564	c======================================================================
	565
	566	do k=1,klev
	567	do i=1,klon
	568	zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
	569	enddo
	570	enddo
	571
	572	c print*,'masse dans ph ',zmasse
[1114]	573	do it=1,nbtr
[541]	574	do k=1,klev
	575	do i=1,klon
	576	d_tr_th(i,k,it)=0.
	577	tr_seri(i,k,it)=max(tr_seri(i,k,it),0.)
	578	tr_seri(i,k,it)=min(tr_seri(i,k,it),1.e10)
	579	enddo
	580	enddo
	581	enddo
	582
[663]	583	if (iflag_thermals.gt.0) then
[644]	584	c print*,'calcul de leffet des thermiques'
[541]	585	nsplit=10
[1114]	586	DO it=1, nbtr
[541]	587	c WRITE(itn,'(i1)') it
	588	c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'conv it='//itn)
	589	c print*,'avant dqthermiquesretro'
	590	c call dump2d(iim,jjm-1,tr_seri(2,1,1),'TR_SERI ')
	591
	592	do isplit=1,nsplit
	593	c Abderr 25 11 02
	594	C Thermiques
[662]	595	c print*,'Avant dans phytrac'
[541]	596	call dqthermcell(klon,klev,pdtphys/nsplit
	597	. ,fm_therm,entr_therm,zmasse
	598	. ,tr_seri(1:klon,1:klev,it),d_tr,ztra_th)
	599
	600	do k=1,klev
	601	do i=1,klon
	602	d_tr(i,k)=pdtphys*d_tr(i,k)/nsplit
	603	d_tr_th(i,k,it)=d_tr_th(i,k,it)+d_tr(i,k)
	604	tr_seri(i,k,it)=max(tr_seri(i,k,it)+d_tr(i,k),0.)
	605	enddo
	606	enddo
	607	enddo ! nsplit
[644]	608	c print*,'apres thermiques'
[541]	609	c call dump2d(iim,jjm-1,d_tr_th(1,1,1),'d_tr_th ')
	610	c do k=1,klev
	611	c print*,'d_tr_th(',k,')=',tr_seri(280,k,1)
	612	c enddo
	613
	614	c WRITE(itn,'(i1)') it
	615	c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'therm it='//itn)
	616	ENDDO ! it
	617	endif ! Thermiques
	618	c print*,'ATTENTION: sdans thermniques'
	619
	620	c======================================================================
[524]	621	c Calcul de l'effet de la couche limite
	622	c======================================================================
	623	c print *,'Avant couchelimite'
[1114]	624	c do it=1,nbtr
[524]	625	c WRITE(itn,'(i1)') it
	626	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
	627	c enddo
	628
	629	if (couchelimite) then
	630
	631	DO k = 1, nlev
	632	DO i = 1, klon
	633	delp(i,k) = paprs(i,k)-paprs(i,k+1)
	634	ENDDO
	635	ENDDO
	636
	637	C maf modif pour tenir compte du cas rnpb + traceur
[1114]	638	DO it=1, nbtr
[959]	639
[1114]	640	IF ( config_inca/='none' ) THEN
	641	IF( pbl_flg(it) == 0 ) CYCLE
	642	END IF
[959]	643
[524]	644	c print *,'it',it,clsol(it)
	645	if (clsol(it)) then ! couche limite avec quantite dans le sol calculee
	646	CALL cltracrn(it, pdtphys, yu1, yv1,
[1067]	647	e cdragh, coefh,t_seri,ftsol,pctsrf,
[524]	648	e tr_seri(1,1,it),trs(1,it),
	649	e paprs, pplay, delp,
	650	e masktr(1,it),fshtr(1,it),hsoltr(it),
	651	e tautr(it),vdeptr(it),
	652	e xlat,
[619]	653	s d_tr_cl(1,1,it),d_trs)
[524]	654	DO k = 1, nlev
	655	DO i = 1, klon
[619]	656	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
[524]	657	ENDDO
	658	ENDDO
	659	c
	660	c Traceur ds sol
	661	c
	662	DO i = 1, klon
	663	trs(i,it) = trs(i,it) + d_trs(i)
	664	END DO
	665	C
	666	C maf provisoire suppression des prints
	667	C WRITE(itn,'(i1)') it
	668	C CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
	669	else ! couche limite avec flux prescrit
	670
[1141]	671	IF (config_inca == 'none') THEN
[524]	672	Cmaf provisoire source / traceur a creer
[1141]	673	DO i=1, klon
	674	source(i,it) = 0.0 ! pas de source, pour l'instant
	675	ENDDO
	676	END IF
	677
[524]	678	CALL cltrac(pdtphys, coefh,t_seri,
[766]	679	s tr_seri(1,1,it), source(:,it),
[524]	680	e paprs, pplay, delp,
[619]	681	s d_tr_cl(1,1,it))
[524]	682	DO k = 1, nlev
	683	DO i = 1, klon
[619]	684	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
[524]	685	ENDDO
	686	ENDDO
	687	Cmaf WRITE(itn,'(i1)') it
	688	cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
	689	endif
	690	ENDDO
	691	c
	692	endif ! couche limite
	693
	694	c print*,'Apres couchelimite'
[1114]	695	c do it=1,nbtr
[524]	696	c WRITE(itn,'(i1)') it
	697	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
	698	c enddo
	699
	700	c======================================================================
	701	c Calcul de l'effet du puits radioactif
	702	c======================================================================
	703
	704	C MAF il faudrait faire une modification pour passer dans radiornpb
	705	c si radio=true mais pour l'instant radiornpb propre au cas rnpb
	706	if(rnpb) then
[793]	707	c print *, 'decroissance radiactive activee'
[524]	708	call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
	709	C
[1114]	710	DO it=1,nbtr
[524]	711	if(radio(it)) then
	712	DO k = 1, nlev
	713	DO i = 1, klon
	714	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
	715	ENDDO
	716	ENDDO
	717	WRITE(itn,'(i1)') it
	718	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn)
	719	endif
	720	ENDDO
	721	c
	722	endif ! rnpb decroissance radioactive
	723	C
	724	c======================================================================
	725	c Calcul de l'effet de la precipitation
	726	c======================================================================
	727
	728	c print*,'LESSIVAGE =',lessivage
	729	IF (lessivage) THEN
	730
	731	c print*,'avant lessivage'
	732
	733	d_tr_lessi_nucl(:,:,:) = 0.
	734	d_tr_lessi_impa(:,:,:) = 0.
	735	flestottr(:,:,:) = 0.
	736	c
	737	c tendance des aerosols nuclees et impactes
	738	c
[1114]	739	DO it = 1, nbtr
[524]	740	IF (aerosol(it)) THEN
	741	DO k = 1, nlev
	742	DO i = 1, klon
	743	d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +
	744	s ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
	745	d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +
	746	s ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
	747	ENDDO
	748	ENDDO
	749	ENDIF
	750	ENDDO
	751	c
	752	c Mises a jour des traceurs + calcul des flux de lessivage
	753	c Mise a jour due a l'impaction et a la nucleation
	754	c
	755	c call dump2d(iim,jjm-1,frac_impa(2:klon-1,10),'FRACIMPA')
	756	c call dump2d(iim,jjm-1,frac_nucl(2:klon-1,10),'FRACNUCL')
	757	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3')
[1114]	758	DO it = 1, nbtr
[524]	759	c print*,'IT=',it,aerosol(it)
	760	IF (aerosol(it)) THEN
	761	c print*,'IT=',it,' On lessive'
	762	DO k = 1, nlev
	763	DO i = 1, klon
	764	tr_seri(i,k,it)=tr_seri(i,k,it)
	765	s frac_impa(i,k)frac_nucl(i,k)
	766	ENDDO
	767	ENDDO
	768	ENDIF
	769	ENDDO
	770	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3B')
	771	c
	772	c Flux lessivage total
	773	c
[1114]	774	DO it = 1, nbtr
[524]	775	DO k = 1, nlev
	776	DO i = 1, klon
	777	flestottr(i,k,it) = flestottr(i,k,it) -
	778	s ( d_tr_lessi_nucl(i,k,it) +
	779	s d_tr_lessi_impa(i,k,it) ) *
	780	s ( paprs(i,k)-paprs(i,k+1) ) /
	781	s (RG * pdtphys)
	782	ENDDO
	783	ENDDO
	784	c
	785	Cmaf WRITE(itn,'(i1)') it
	786	Cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
	787	ENDDO
	788	c
	789	c print*,'apres lessivage'
	790	ENDIF
	791	Cc
	792	DO k = 1, nlev
	793	DO i = 1, klon
	794	fluxrn(i,k) = flestottr(i,k,1)
	795	fluxpb(i,k) = flestottr(i,k,2)
	796	rn(i,k) = tr_seri(i,k,1)
	797	pb(i,k) = tr_seri(i,k,2)
	798	pbnucl(i,k)=d_tr_lessi_nucl(i,k,2)
	799	pbimpa(i,k)=d_tr_lessi_impa(i,k,2)
	800	ENDDO
	801	ENDDO
	802
	803	c=============================================================
	804	c Ecriture des sorties
	805	c=============================================================
	806
	807	#ifdef CPP_IOIPSL
	808	#include "write_histrac.h"
	809	#endif
	810
	811	c=============================================================
	812
	813	if (lafin) then
	814	print*, 'c est la fin de la physique'
[776]	815	call Gather(trs(:,1),trs_tmp)
[766]	816	c$OMP MASTER
[776]	817	if (is_mpi_root) then
[766]	818
	819	open (99,file='restarttrac', form='formatted')
[776]	820	do i=1,klon_glo
[766]	821	write(99,*) trs_tmp(i)
	822	enddo
	823	PRINT*, 'Ecriture du fichier restarttrac'
	824	close(99)
	825	endif
	826	c$OMP END MASTER
[524]	827	else
[644]	828	c print*, 'physique pas fini'
[524]	829	endif
	830
	831
	832	RETURN
	833	END

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

Download in other formats:

Original Format