Context Navigation

phytrac.F @ 1085

Last change on this file since 1085 was 1067, checked in by Laurent Fairhead, 16 years ago

Modifications lie au premier niveau du modele pour la diffusion turbulent

du vent.

Preparation pour un couplage des courrant oceaniques.

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

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

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

Download in other formats:

Original Format