Context Navigation

phytrac.F @ 1219

Last change on this file since 1219 was 1146, checked in by Laurent Fairhead, 16 years ago
Réintegration dans le tronc des modifications issues de la branche LMDZ-dev comprises entre la révision 1074 et 1145 Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes résultats entre la trunk et la dev LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 26.5 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,
[959]	59	I tau_inca,
	60	I piz_inca,
	61	I cg_inca,
	62	I ccm,
	63	I rfname,
[524]	64	O tr_seri)
	65
	66	USE ioipsl
[766]	67	USE dimphy
[1146]	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)
[1146]	138	CHARACTER(len=8) :: solsym(nbtr)
[864]	139	integer la
	140	REAL :: tau_inca(klon,klev,9,2)
	141	REAL :: piz_inca(klon,klev,9,2)
	142	REAL :: cg_inca(klon,klev,9,2)
	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
[1146]	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
[1146]	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,
[864]	460	$ tau_inca,
	461	$ piz_inca,
	462	$ cg_inca,
	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	$ .false., !wrhstts
	500	$ hbuf, !hbuf
	501	$ obuf, !obuf
	502	$ iip1, !nx
[766]	503	$ jjp1, !ny
	504	$ source,
	505	$ solsym)
[682]	506
	507
[766]	508	call VTe(VTinca)
	509	call VTb(VTphysiq)
[524]	510	#endif
[959]	511	END IF ! config_inca
[524]	512	c======================================================================
	513	c Calcul de l'effet de la convection
	514	c======================================================================
	515	c print*,'Avant convection'
[1146]	516	do it=1,nbtr
[928]	517	WRITE(itn,'(i2)') it
[524]	518	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
[928]	519	enddo
[524]	520
	521	if (convection) then
	522
	523	c print*,'Pas de temps dans phytrac : ',pdtphys
[1146]	524	DO it=1, nbtr
[959]	525
[1146]	526	IF ( config_inca/='none') THEN
	527	IF ( conv_flg(it) == 0 ) CYCLE
	528	END IF
[959]	529
[928]	530	if (iflag_con.lt.2) then
	531	d_tr_cv=0.
	532	else if (iflag_con.eq.2) then
[619]	533	c tiedke
[524]	534	CALL nflxtr(pdtphys, pmfu, pmfd, pen_u, pde_u, pen_d, pde_d,
[619]	535	. pplay, paprs, tr_seri(1,1,it), d_tr_cv(1,1,it))
[928]	536	else
[619]	537	c KE
	538	call cvltr(pdtphys, da, phi, mp, paprs,pplay, tr_seri(1,1,it),
	539	. upwd,dnwd,d_tr_cv(1,1,it))
	540	endif
	541
[928]	542	DO k = 1, nlev
	543	DO i = 1, klon
[619]	544	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cv(i,k,it)
[928]	545	ENDDO
	546	ENDDO
[959]	547
	548	IF (config_inca == 'none') THEN
	549	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '//itn)
	550	ELSE
	551	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'convection it = '
	552	. //solsym(it))
	553	END IF
	554
[524]	555	ENDDO
	556
	557	endif ! convection
	558	c print*,'Apres convection'
[1146]	559	c do it=1,nbtr
[524]	560	c WRITE(itn,'(i1)') it
	561	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant conv'//itn)
	562	c enddo
	563
[541]	564
[524]	565	c======================================================================
[541]	566	c Calcul de l'effet des thermiques
	567	c======================================================================
	568
	569	do k=1,klev
	570	do i=1,klon
	571	zmasse(i,k)=(paprs(i,k)-paprs(i,k+1))/rg
	572	enddo
	573	enddo
	574
	575	c print*,'masse dans ph ',zmasse
[1146]	576	do it=1,nbtr
[541]	577	do k=1,klev
	578	do i=1,klon
	579	d_tr_th(i,k,it)=0.
	580	tr_seri(i,k,it)=max(tr_seri(i,k,it),0.)
	581	tr_seri(i,k,it)=min(tr_seri(i,k,it),1.e10)
	582	enddo
	583	enddo
	584	enddo
	585
[663]	586	if (iflag_thermals.gt.0) then
[644]	587	c print*,'calcul de leffet des thermiques'
[541]	588	nsplit=10
[1146]	589	DO it=1, nbtr
[541]	590	c WRITE(itn,'(i1)') it
	591	c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'conv it='//itn)
	592	c print*,'avant dqthermiquesretro'
	593	c call dump2d(iim,jjm-1,tr_seri(2,1,1),'TR_SERI ')
	594
	595	do isplit=1,nsplit
	596	c Abderr 25 11 02
	597	C Thermiques
[662]	598	c print*,'Avant dans phytrac'
[541]	599	call dqthermcell(klon,klev,pdtphys/nsplit
	600	. ,fm_therm,entr_therm,zmasse
	601	. ,tr_seri(1:klon,1:klev,it),d_tr,ztra_th)
	602
	603	do k=1,klev
	604	do i=1,klon
	605	d_tr(i,k)=pdtphys*d_tr(i,k)/nsplit
	606	d_tr_th(i,k,it)=d_tr_th(i,k,it)+d_tr(i,k)
	607	tr_seri(i,k,it)=max(tr_seri(i,k,it)+d_tr(i,k),0.)
	608	enddo
	609	enddo
	610	enddo ! nsplit
[644]	611	c print*,'apres thermiques'
[541]	612	c call dump2d(iim,jjm-1,d_tr_th(1,1,1),'d_tr_th ')
	613	c do k=1,klev
	614	c print*,'d_tr_th(',k,')=',tr_seri(280,k,1)
	615	c enddo
	616
	617	c WRITE(itn,'(i1)') it
	618	c CALL minmaxqfi(tr_seri(1,1,it),1.e10,-1.e33,'therm it='//itn)
	619	ENDDO ! it
	620	endif ! Thermiques
	621	c print*,'ATTENTION: sdans thermniques'
	622
	623	c======================================================================
[524]	624	c Calcul de l'effet de la couche limite
	625	c======================================================================
	626	c print *,'Avant couchelimite'
[1146]	627	c do it=1,nbtr
[524]	628	c WRITE(itn,'(i1)') it
	629	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
	630	c enddo
	631
	632	if (couchelimite) then
	633
	634	DO k = 1, nlev
	635	DO i = 1, klon
	636	delp(i,k) = paprs(i,k)-paprs(i,k+1)
	637	ENDDO
	638	ENDDO
	639
	640	C maf modif pour tenir compte du cas rnpb + traceur
[1146]	641	DO it=1, nbtr
[959]	642
[1146]	643	IF ( config_inca/='none' ) THEN
	644	IF( pbl_flg(it) == 0 ) CYCLE
	645	END IF
[959]	646
[524]	647	c print *,'it',it,clsol(it)
	648	if (clsol(it)) then ! couche limite avec quantite dans le sol calculee
	649	CALL cltracrn(it, pdtphys, yu1, yv1,
[1067]	650	e cdragh, coefh,t_seri,ftsol,pctsrf,
[524]	651	e tr_seri(1,1,it),trs(1,it),
	652	e paprs, pplay, delp,
	653	e masktr(1,it),fshtr(1,it),hsoltr(it),
	654	e tautr(it),vdeptr(it),
	655	e xlat,
[619]	656	s d_tr_cl(1,1,it),d_trs)
[524]	657	DO k = 1, nlev
	658	DO i = 1, klon
[619]	659	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
[524]	660	ENDDO
	661	ENDDO
	662	c
	663	c Traceur ds sol
	664	c
	665	DO i = 1, klon
	666	trs(i,it) = trs(i,it) + d_trs(i)
	667	END DO
	668	C
	669	C maf provisoire suppression des prints
	670	C WRITE(itn,'(i1)') it
	671	C CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracrn it='//itn)
	672	else ! couche limite avec flux prescrit
	673
[1146]	674	IF (config_inca == 'none') THEN
[524]	675	Cmaf provisoire source / traceur a creer
[1146]	676	DO i=1, klon
	677	source(i,it) = 0.0 ! pas de source, pour l'instant
	678	ENDDO
	679	END IF
	680
[524]	681	CALL cltrac(pdtphys, coefh,t_seri,
[766]	682	s tr_seri(1,1,it), source(:,it),
[524]	683	e paprs, pplay, delp,
[619]	684	s d_tr_cl(1,1,it))
[524]	685	DO k = 1, nlev
	686	DO i = 1, klon
[619]	687	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_cl(i,k,it)
[524]	688	ENDDO
	689	ENDDO
	690	Cmaf WRITE(itn,'(i1)') it
	691	cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'cltracn it='//itn)
	692	endif
	693	ENDDO
	694	c
	695	endif ! couche limite
	696
	697	c print*,'Apres couchelimite'
[1146]	698	c do it=1,nbtr
[524]	699	c WRITE(itn,'(i1)') it
	700	c call diagtracphy(tr_seri(:,:,it),paprs,'Avant CL '//itn)
	701	c enddo
	702
	703	c======================================================================
	704	c Calcul de l'effet du puits radioactif
	705	c======================================================================
	706
	707	C MAF il faudrait faire une modification pour passer dans radiornpb
	708	c si radio=true mais pour l'instant radiornpb propre au cas rnpb
	709	if(rnpb) then
[793]	710	c print *, 'decroissance radiactive activee'
[524]	711	call radiornpb (tr_seri,pdtphys,tautr,d_tr_dec)
	712	C
[1146]	713	DO it=1,nbtr
[524]	714	if(radio(it)) then
	715	DO k = 1, nlev
	716	DO i = 1, klon
	717	tr_seri(i,k,it) = tr_seri(i,k,it) + d_tr_dec(i,k,it)
	718	ENDDO
	719	ENDDO
	720	WRITE(itn,'(i1)') it
	721	CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'puits rn it='//itn)
	722	endif
	723	ENDDO
	724	c
	725	endif ! rnpb decroissance radioactive
	726	C
	727	c======================================================================
	728	c Calcul de l'effet de la precipitation
	729	c======================================================================
	730
	731	c print*,'LESSIVAGE =',lessivage
	732	IF (lessivage) THEN
	733
	734	c print*,'avant lessivage'
	735
	736	d_tr_lessi_nucl(:,:,:) = 0.
	737	d_tr_lessi_impa(:,:,:) = 0.
	738	flestottr(:,:,:) = 0.
	739	c
	740	c tendance des aerosols nuclees et impactes
	741	c
[1146]	742	DO it = 1, nbtr
[524]	743	IF (aerosol(it)) THEN
	744	DO k = 1, nlev
	745	DO i = 1, klon
	746	d_tr_lessi_nucl(i,k,it) = d_tr_lessi_nucl(i,k,it) +
	747	s ( 1 - frac_nucl(i,k) )*tr_seri(i,k,it)
	748	d_tr_lessi_impa(i,k,it) = d_tr_lessi_impa(i,k,it) +
	749	s ( 1 - frac_impa(i,k) )*tr_seri(i,k,it)
	750	ENDDO
	751	ENDDO
	752	ENDIF
	753	ENDDO
	754	c
	755	c Mises a jour des traceurs + calcul des flux de lessivage
	756	c Mise a jour due a l'impaction et a la nucleation
	757	c
	758	c call dump2d(iim,jjm-1,frac_impa(2:klon-1,10),'FRACIMPA')
	759	c call dump2d(iim,jjm-1,frac_nucl(2:klon-1,10),'FRACNUCL')
	760	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3')
[1146]	761	DO it = 1, nbtr
[524]	762	c print*,'IT=',it,aerosol(it)
	763	IF (aerosol(it)) THEN
	764	c print*,'IT=',it,' On lessive'
	765	DO k = 1, nlev
	766	DO i = 1, klon
	767	tr_seri(i,k,it)=tr_seri(i,k,it)
	768	s frac_impa(i,k)frac_nucl(i,k)
	769	ENDDO
	770	ENDDO
	771	ENDIF
	772	ENDDO
	773	c call dump2d(iim,jjm-1,tr_seri(2:klon-1,10,3),'TRACEUR3B')
	774	c
	775	c Flux lessivage total
	776	c
[1146]	777	DO it = 1, nbtr
[524]	778	DO k = 1, nlev
	779	DO i = 1, klon
	780	flestottr(i,k,it) = flestottr(i,k,it) -
	781	s ( d_tr_lessi_nucl(i,k,it) +
	782	s d_tr_lessi_impa(i,k,it) ) *
	783	s ( paprs(i,k)-paprs(i,k+1) ) /
	784	s (RG * pdtphys)
	785	ENDDO
	786	ENDDO
	787	c
	788	Cmaf WRITE(itn,'(i1)') it
	789	Cmaf CALL minmaxqfi(tr_seri(1,1,it),0.,1.e33,'tr(lessi) it='//itn)
	790	ENDDO
	791	c
	792	c print*,'apres lessivage'
	793	ENDIF
	794	Cc
	795	DO k = 1, nlev
	796	DO i = 1, klon
	797	fluxrn(i,k) = flestottr(i,k,1)
	798	fluxpb(i,k) = flestottr(i,k,2)
	799	rn(i,k) = tr_seri(i,k,1)
	800	pb(i,k) = tr_seri(i,k,2)
	801	pbnucl(i,k)=d_tr_lessi_nucl(i,k,2)
	802	pbimpa(i,k)=d_tr_lessi_impa(i,k,2)
	803	ENDDO
	804	ENDDO
	805
	806	c=============================================================
	807	c Ecriture des sorties
	808	c=============================================================
	809
	810	#ifdef CPP_IOIPSL
	811	#include "write_histrac.h"
	812	#endif
	813
	814	c=============================================================
	815
	816	if (lafin) then
	817	print*, 'c est la fin de la physique'
[776]	818	call Gather(trs(:,1),trs_tmp)
[766]	819	c$OMP MASTER
[776]	820	if (is_mpi_root) then
[766]	821
	822	open (99,file='restarttrac', form='formatted')
[776]	823	do i=1,klon_glo
[766]	824	write(99,*) trs_tmp(i)
	825	enddo
	826	PRINT*, 'Ecriture du fichier restarttrac'
	827	close(99)
	828	endif
	829	c$OMP END MASTER
[524]	830	else
[644]	831	c print*, 'physique pas fini'
[524]	832	endif
	833
	834
	835	RETURN
	836	END

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

Download in other formats:

Original Format