Context Navigation

bilan_dyn_p.F @ 1368

Last change on this file since 1368 was 1279, checked in by Laurent Fairhead, 15 years ago
Merged LMDZ4-dev branch changes r1241:1278 into the trunk Running trunk and LMDZ4-dev in LMDZOR configuration on local machine (sequential) and SX8 (4-proc) yields identical results (restart and restartphy are identical binarily) Log history from r1241 to r1278 is available by switching to source:LMDZ4/branches/LMDZ4-dev-20091210
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 19.9 KB

Rev	Line
[630]	1	!
[1279]	2	! $Id: bilan_dyn_p.F 1279 2009-12-10 09:02:56Z idelkadi $
[630]	3	!
	4	SUBROUTINE bilan_dyn_p (ntrac,dt_app,dt_cum,
	5	s ps,masse,pk,flux_u,flux_v,teta,phi,ucov,vcov,trac)
	6
	7	c AFAIRE
	8	c Prevoir en champ nq+1 le diagnostique de l'energie
	9	c en faisant Qzon=Cv T + L * ...
	10	c vQ..A=Cp T + L * ...
	11
[1279]	12	#ifdef CPP_IOIPSL
[630]	13	USE IOIPSL
[1279]	14	#endif
[630]	15	USE parallel
	16	USE mod_hallo
	17	use misc_mod
	18	use write_field
	19	IMPLICIT NONE
	20
	21	#include "dimensions.h"
	22	#include "paramet.h"
	23	#include "comconst.h"
	24	#include "comvert.h"
	25	#include "comgeom2.h"
	26	#include "temps.h"
	27	#include "iniprint.h"
	28
	29	c====================================================================
	30	c
	31	c Sous-programme consacre à des diagnostics dynamiques de base
	32	c
	33	c
	34	c De facon generale, les moyennes des scalaires Q sont ponderees par
	35	c la masse.
	36	c
	37	c Les flux de masse sont eux simplement moyennes.
	38	c
	39	c====================================================================
	40
	41	c Arguments :
	42	c ===========
	43
	44	integer ntrac
	45	real dt_app,dt_cum
	46	real ps(iip1,jjp1)
	47	real masse(iip1,jjp1,llm),pk(iip1,jjp1,llm)
	48	real flux_u(iip1,jjp1,llm)
	49	real flux_v(iip1,jjm,llm)
	50	real teta(iip1,jjp1,llm)
	51	real phi(iip1,jjp1,llm)
	52	real ucov(iip1,jjp1,llm)
	53	real vcov(iip1,jjm,llm)
	54	real trac(iip1,jjp1,llm,ntrac)
	55
	56	c Local :
	57	c =======
	58
	59	integer icum,ncum
	60	logical first
	61	real zz,zqy,zfactv(jjm,llm)
	62
	63	integer nQ
	64	parameter (nQ=7)
	65
	66
	67	cym character*6 nom(nQ)
	68	cym character*6 unites(nQ)
	69	character*6,save :: nom(nQ)
	70	character*6,save :: unites(nQ)
	71
	72	character*10 file
	73	integer ifile
	74	parameter (ifile=4)
	75
	76	integer itemp,igeop,iecin,iang,iu,iovap,iun
	77	integer i_sortie
	78
	79	save first,icum,ncum
	80	save itemp,igeop,iecin,iang,iu,iovap,iun
	81	save i_sortie
	82
	83	real time
	84	integer itau
	85	save time,itau
	86	data time,itau/0.,0/
	87
	88	data first/.true./
	89	data itemp,igeop,iecin,iang,iu,iovap,iun/1,2,3,4,5,6,7/
	90	data i_sortie/1/
	91
	92	real ww
	93
	94	c variables dynamiques intermédiaires
	95	REAL vcont(iip1,jjm,llm),ucont(iip1,jjp1,llm)
	96	REAL ang(iip1,jjp1,llm),unat(iip1,jjp1,llm)
	97	REAL massebx(iip1,jjp1,llm),masseby(iip1,jjm,llm)
	98	REAL vorpot(iip1,jjm,llm)
	99	REAL w(iip1,jjp1,llm),ecin(iip1,jjp1,llm),convm(iip1,jjp1,llm)
	100	REAL bern(iip1,jjp1,llm)
	101
	102	c champ contenant les scalaires advectés.
	103	real Q(iip1,jjp1,llm,nQ)
	104
	105	c champs cumulés
	106	real ps_cum(iip1,jjp1)
	107	real masse_cum(iip1,jjp1,llm)
	108	real flux_u_cum(iip1,jjp1,llm)
	109	real flux_v_cum(iip1,jjm,llm)
	110	real Q_cum(iip1,jjp1,llm,nQ)
	111	real flux_uQ_cum(iip1,jjp1,llm,nQ)
	112	real flux_vQ_cum(iip1,jjm,llm,nQ)
	113	real flux_wQ_cum(iip1,jjp1,llm,nQ)
	114	real dQ(iip1,jjp1,llm,nQ)
	115
	116	save ps_cum,masse_cum,flux_u_cum,flux_v_cum
	117	save Q_cum,flux_uQ_cum,flux_vQ_cum
	118
	119	c champs de tansport en moyenne zonale
	120	integer ntr,itr
	121	parameter (ntr=5)
	122
	123	cym character*10 znom(ntr,nQ)
	124	cym character*20 znoml(ntr,nQ)
	125	cym character*10 zunites(ntr,nQ)
	126	character*10,save :: znom(ntr,nQ)
	127	character*20,save :: znoml(ntr,nQ)
	128	character*10,save :: zunites(ntr,nQ)
	129
	130	integer iave,itot,immc,itrs,istn
	131	data iave,itot,immc,itrs,istn/1,2,3,4,5/
	132	character*3 ctrs(ntr)
	133	data ctrs/' ','TOT','MMC','TRS','STN'/
	134
	135	real zvQ(jjm,llm,ntr,nQ),zvQtmp(jjm,llm)
	136	real zavQ(jjm,ntr,nQ),psiQ(jjm,llm+1,nQ)
	137	real zmasse(jjm,llm),zamasse(jjm)
	138
	139	real zv(jjm,llm),psi(jjm,llm+1)
	140
	141	integer i,j,l,iQ
	142
	143
	144	c Initialisation du fichier contenant les moyennes zonales.
	145	c ---------------------------------------------------------
	146
	147	character*10 infile
	148
	149	integer fileid
	150	integer thoriid, zvertiid
	151	save fileid
	152
	153	integer ndex3d(jjm*llm)
	154
	155	C Variables locales
	156	C
	157	integer tau0
	158	real zjulian
	159	character*3 str
	160	character*10 ctrac
	161	integer ii,jj
	162	integer zan, dayref
	163	C
	164	real rlong(jjm),rlatg(jjm)
	165	integer :: jjb,jje,jjn,ijb,ije
	166	type(Request) :: Req
	167
[985]	168	! definition du domaine d'ecriture pour le rebuild
[630]	169
[985]	170	INTEGER,DIMENSION(1) :: ddid
	171	INTEGER,DIMENSION(1) :: dsg
	172	INTEGER,DIMENSION(1) :: dsl
	173	INTEGER,DIMENSION(1) :: dpf
	174	INTEGER,DIMENSION(1) :: dpl
	175	INTEGER,DIMENSION(1) :: dhs
	176	INTEGER,DIMENSION(1) :: dhe
	177
	178	INTEGER :: bilan_dyn_domain_id
	179
	180
[630]	181	c=====================================================================
	182	c Initialisation
	183	c=====================================================================
[764]	184	ndex3d=0
[630]	185	if (adjust) return
	186
	187	time=time+dt_app
	188	itau=itau+1
	189
	190	if (first) then
	191
	192
	193	icum=0
	194	c initialisation des fichiers
	195	first=.false.
	196	c ncum est la frequence de stokage en pas de temps
	197	ncum=dt_cum/dt_app
	198	if (abs(ncumdt_app-dt_cum).gt.1.e-5dt_app) then
	199	WRITE(lunout,*)
	200	. 'Pb : le pas de cumule doit etre multiple du pas'
	201	WRITE(lunout,*)'dt_app=',dt_app
	202	WRITE(lunout,*)'dt_cum=',dt_cum
	203	stop
	204	endif
	205
	206	if (i_sortie.eq.1) then
	207	file='dynzon'
	208	if (mpi_rank==0) then
	209	call inigrads(ifile,1
	210	s ,0.,180./pi,0.,0.,jjm,rlatv,-90.,90.,180./pi
	211	s ,llm,presnivs,1.
	212	s ,dt_cum,file,'dyn_zon ')
	213	endif
	214	endif
	215
	216	nom(itemp)='T'
	217	nom(igeop)='gz'
	218	nom(iecin)='K'
	219	nom(iang)='ang'
	220	nom(iu)='u'
	221	nom(iovap)='ovap'
	222	nom(iun)='un'
	223
	224	unites(itemp)='K'
	225	unites(igeop)='m2/s2'
	226	unites(iecin)='m2/s2'
	227	unites(iang)='ang'
	228	unites(iu)='m/s'
	229	unites(iovap)='kg/kg'
	230	unites(iun)='un'
	231
	232
	233	c Initialisation du fichier contenant les moyennes zonales.
	234	c ---------------------------------------------------------
	235
	236	infile='dynzon'
	237
	238	zan = annee_ref
	239	dayref = day_ref
	240	CALL ymds2ju(zan, 1, dayref, 0.0, zjulian)
	241	tau0 = itau_dyn
	242
	243	rlong=0.
	244	rlatg=rlatv*180./pi
	245
	246	jjb=jj_begin
	247	jje=jj_end
	248	jjn=jj_nb
[985]	249	IF (pole_sud) THEN
	250	jjn=jj_nb-1
	251	jje=jj_end-1
	252	ENDIF
	253
	254	ddid=(/ 2 /)
	255	dsg=(/ jjm /)
	256	dsl=(/ jjn /)
	257	dpf=(/ jjb /)
	258	dpl=(/ jje /)
	259	dhs=(/ 0 /)
	260	dhe=(/ 0 /)
	261
	262	call flio_dom_set(mpi_size,mpi_rank,ddid,dsg,dsl,dpf,dpl,dhs,dhe,
	263	. 'box',bilan_dyn_domain_id)
[630]	264
[985]	265	call histbeg(trim(infile),
[630]	266	. 1, rlong(jjb:jje), jjn, rlatg(jjb:jje),
	267	. 1, 1, 1, jjn,
[985]	268	. tau0, zjulian, dt_cum, thoriid, fileid,
	269	. bilan_dyn_domain_id)
[630]	270
	271	C
	272	C Appel a histvert pour la grille verticale
	273	C
	274	call histvert(fileid, 'presnivs', 'Niveaux sigma','mb',
	275	. llm, presnivs, zvertiid)
	276	C
	277	C Appels a histdef pour la definition des variables a sauvegarder
	278	do iQ=1,nQ
	279	do itr=1,ntr
	280	if(itr.eq.1) then
	281	znom(itr,iQ)=nom(iQ)
	282	znoml(itr,iQ)=nom(iQ)
	283	zunites(itr,iQ)=unites(iQ)
	284	else
	285	znom(itr,iQ)=ctrs(itr)//'v'//nom(iQ)
	286	znoml(itr,iQ)='transport : v * '//nom(iQ)//' '//ctrs(itr)
	287	zunites(itr,iQ)='m/s * '//unites(iQ)
	288	endif
	289	enddo
	290	enddo
	291
	292	c Declarations des champs avec dimension verticale
	293	c print*,'1HISTDEF'
	294	do iQ=1,nQ
	295	do itr=1,ntr
	296	IF (prt_level > 5)
	297	. WRITE(lunout,*)'var ',itr,iQ
	298	. ,znom(itr,iQ),znoml(itr,iQ),zunites(itr,iQ)
	299	call histdef(fileid,znom(itr,iQ),znoml(itr,iQ),
	300	. zunites(itr,iQ),1,jjn,thoriid,llm,1,llm,zvertiid,
	301	. 32,'ave(X)',dt_cum,dt_cum)
	302	enddo
	303	c Declarations pour les fonctions de courant
	304	c print*,'2HISTDEF'
	305	call histdef(fileid,'psi'//nom(iQ)
	306	. ,'stream fn. '//znoml(itot,iQ),
	307	. zunites(itot,iQ),1,jjn,thoriid,llm,1,llm,zvertiid,
	308	. 32,'ave(X)',dt_cum,dt_cum)
	309	enddo
	310
	311
	312	c Declarations pour les champs de transport d'air
	313	c print*,'3HISTDEF'
	314	call histdef(fileid, 'masse', 'masse',
	315	. 'kg', 1, jjn, thoriid, llm, 1, llm, zvertiid,
	316	. 32, 'ave(X)', dt_cum, dt_cum)
	317	call histdef(fileid, 'v', 'v',
	318	. 'm/s', 1, jjn, thoriid, llm, 1, llm, zvertiid,
	319	. 32, 'ave(X)', dt_cum, dt_cum)
	320	c Declarations pour les fonctions de courant
	321	c print*,'4HISTDEF'
	322	call histdef(fileid,'psi','stream fn. MMC ','mega t/s',
	323	. 1,jjn,thoriid,llm,1,llm,zvertiid,
	324	. 32,'ave(X)',dt_cum,dt_cum)
	325
	326
	327	c Declaration des champs 1D de transport en latitude
	328	c print*,'5HISTDEF'
	329	do iQ=1,nQ
	330	do itr=2,ntr
	331	call histdef(fileid,'a'//znom(itr,iQ),znoml(itr,iQ),
	332	. zunites(itr,iQ),1,jjn,thoriid,1,1,1,-99,
	333	. 32,'ave(X)',dt_cum,dt_cum)
	334	enddo
	335	enddo
	336
	337
	338	c print*,'8HISTDEF'
	339	CALL histend(fileid)
	340
	341
	342	endif
	343
	344
	345	c=====================================================================
	346	c Calcul des champs dynamiques
	347	c ----------------------------
	348
	349	jjb=jj_begin
	350	jje=jj_end
	351
	352	c énergie cinétique
	353	ucont(:,jjb:jje,:)=0
	354
	355	call Register_Hallo(ucov,ip1jmp1,llm,1,1,1,1,Req)
	356	call Register_Hallo(vcov,ip1jm,llm,1,1,1,1,Req)
	357	call SendRequest(Req)
	358	call WaitRequest(Req)
	359
	360	CALL covcont_p(llm,ucov,vcov,ucont,vcont)
	361	CALL enercin_p(vcov,ucov,vcont,ucont,ecin)
	362
	363	c moment cinétique
	364	do l=1,llm
	365	ang(:,jjb:jje,l)=ucov(:,jjb:jje,l)+constang(:,jjb:jje)
	366	unat(:,jjb:jje,l)=ucont(:,jjb:jje,l)*cu(:,jjb:jje)
	367	enddo
	368
	369	Q(:,jjb:jje,:,itemp)=teta(:,jjb:jje,:)*pk(:,jjb:jje,:)/cpp
	370	Q(:,jjb:jje,:,igeop)=phi(:,jjb:jje,:)
	371	Q(:,jjb:jje,:,iecin)=ecin(:,jjb:jje,:)
	372	Q(:,jjb:jje,:,iang)=ang(:,jjb:jje,:)
	373	Q(:,jjb:jje,:,iu)=unat(:,jjb:jje,:)
[764]	374	Q(:,jjb:jje,:,iovap)=trac(:,jjb:jje,:,1)
[630]	375	Q(:,jjb:jje,:,iun)=1.
	376
	377
	378	c=====================================================================
	379	c Cumul
	380	c=====================================================================
	381	c
	382	if(icum.EQ.0) then
[985]	383	jjb=jj_begin
	384	jje=jj_end
	385
	386	ps_cum(:,jjb:jje)=0.
	387	masse_cum(:,jjb:jje,:)=0.
	388	flux_u_cum(:,jjb:jje,:)=0.
	389	Q_cum(:,jjb:jje,:,:)=0.
	390	flux_uQ_cum(:,jjb:jje,:,:)=0.
	391	flux_v_cum(:,jjb:jje,:)=0.
	392	if (pole_sud) jje=jj_end-1
	393	flux_v_cum(:,jjb:jje,:)=0.
	394	flux_vQ_cum(:,jjb:jje,:,:)=0.
[630]	395	endif
	396
	397	IF (prt_level > 5)
	398	. WRITE(lunout,*)'dans bilan_dyn ',icum,'->',icum+1
	399	icum=icum+1
	400
	401	c accumulation des flux de masse horizontaux
[985]	402	jjb=jj_begin
	403	jje=jj_end
	404
	405	ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)+ps(:,jjb:jje)
	406	masse_cum(:,jjb:jje,:)=masse_cum(:,jjb:jje,:)+masse(:,jjb:jje,:)
	407	flux_u_cum(:,jjb:jje,:)=flux_u_cum(:,jjb:jje,:)
	408	. +flux_u(:,jjb:jje,:)
	409	if (pole_sud) jje=jj_end-1
	410	flux_v_cum(:,jjb:jje,:)=flux_v_cum(:,jjb:jje,:)
	411	. +flux_v(:,jjb:jje,:)
	412
	413	jjb=jj_begin
	414	jje=jj_end
	415
[630]	416	do iQ=1,nQ
	417	Q_cum(:,jjb:jje,:,iQ)=Q_cum(:,jjb:jje,:,iQ)
	418	. +Q(:,jjb:jje,:,iQ)*masse(:,jjb:jje,:)
	419	enddo
	420
	421	c=====================================================================
	422	c FLUX ET TENDANCES
	423	c=====================================================================
	424
	425	c Flux longitudinal
	426	c -----------------
	427	do iQ=1,nQ
	428	do l=1,llm
	429	do j=jjb,jje
	430	do i=1,iim
	431	flux_uQ_cum(i,j,l,iQ)=flux_uQ_cum(i,j,l,iQ)
	432	s +flux_u(i,j,l)0.5(Q(i,j,l,iQ)+Q(i+1,j,l,iQ))
	433	enddo
	434	flux_uQ_cum(iip1,j,l,iQ)=flux_uQ_cum(1,j,l,iQ)
	435	enddo
	436	enddo
	437	enddo
	438
	439	c flux méridien
	440	c -------------
	441	do iQ=1,nQ
[985]	442	call Register_Hallo(Q(1,1,1,iQ),ip1jmp1,llm,0,1,1,0,Req)
	443	enddo
	444	call SendRequest(Req)
	445	call WaitRequest(Req)
	446
	447	jjb=jj_begin
	448	jje=jj_end
	449	if (pole_sud) jje=jj_end-1
	450
	451	do iQ=1,nQ
[630]	452	do l=1,llm
	453	do j=jjb,jje
	454	do i=1,iip1
	455	flux_vQ_cum(i,j,l,iQ)=flux_vQ_cum(i,j,l,iQ)
	456	s +flux_v(i,j,l)0.5(Q(i,j,l,iQ)+Q(i,j+1,l,iQ))
	457	enddo
	458	enddo
	459	enddo
	460	enddo
	461
	462
	463	c tendances
	464	c ---------
	465
	466	c convergence horizontale
	467	call Register_Hallo(flux_uQ_cum,ip1jmp1,llm,2,2,2,2,Req)
	468	call Register_Hallo(flux_vQ_cum,ip1jm,llm,2,2,2,2,Req)
	469	call SendRequest(Req)
	470	call WaitRequest(Req)
	471
	472	call convflu_p(flux_uQ_cum,flux_vQ_cum,llm*nQ,dQ)
	473
	474	c calcul de la vitesse verticale
	475	call Register_Hallo(flux_u_cum,ip1jmp1,llm,2,2,2,2,Req)
	476	call Register_Hallo(flux_v_cum,ip1jm,llm,2,2,2,2,Req)
	477	call SendRequest(Req)
	478	call WaitRequest(Req)
	479
	480	call convmas_p(flux_u_cum,flux_v_cum,convm)
	481	CALL vitvert_p(convm,w)
	482
[985]	483	jjb=jj_begin
	484	jje=jj_end
	485
[630]	486	do iQ=1,nQ
	487	do l=1,llm-1
	488	do j=jjb,jje
	489	do i=1,iip1
	490	ww=-0.5w(i,j,l+1)(Q(i,j,l,iQ)+Q(i,j,l+1,iQ))
	491	dQ(i,j,l ,iQ)=dQ(i,j,l ,iQ)-ww
	492	dQ(i,j,l+1,iQ)=dQ(i,j,l+1,iQ)+ww
	493	enddo
	494	enddo
	495	enddo
	496	enddo
	497	IF (prt_level > 5)
	498	. WRITE(lunout,*)'Apres les calculs fait a chaque pas'
	499	c=====================================================================
	500	c PAS DE TEMPS D'ECRITURE
	501	c=====================================================================
	502	if (icum.eq.ncum) then
	503	c=====================================================================
	504
	505	IF (prt_level > 5)
	506	. WRITE(lunout,*)'Pas d ecriture'
	507
	508	c Normalisation
	509	do iQ=1,nQ
	510	Q_cum(:,jjb:jje,:,iQ)=Q_cum(:,jjb:jje,:,iQ)
	511	. /masse_cum(:,jjb:jje,:)
	512	enddo
	513	zz=1./float(ncum)
	514
[985]	515	jjb=jj_begin
	516	jje=jj_end
[630]	517
[985]	518	ps_cum(:,jjb:jje)=ps_cum(:,jjb:jje)*zz
	519	masse_cum(:,jjb:jje,:)=masse_cum(:,jjb:jje,:)*zz
	520	flux_u_cum(:,jjb:jje,:)=flux_u_cum(:,jjb:jje,:)*zz
	521	flux_uQ_cum(:,jjb:jje,:,:)=flux_uQ_cum(:,jjb:jje,:,:)*zz
	522	dQ(:,jjb:jje,:,:)=dQ(:,jjb:jje,:,:)*zz
	523
	524	IF (pole_sud) jje=jj_end-1
	525	flux_v_cum(:,jjb:jje,:)=flux_v_cum(:,jjb:jje,:)*zz
	526	flux_vQ_cum(:,jjb:jje,:,:)=flux_vQ_cum(:,jjb:jje,:,:)*zz
	527
	528	jjb=jj_begin
	529	jje=jj_end
	530
	531
[630]	532	c A retravailler eventuellement
	533	c division de dQ par la masse pour revenir aux bonnes grandeurs
	534	do iQ=1,nQ
	535	dQ(:,jjb:jje,:,iQ)=dQ(:,jjb:jje,:,iQ)/masse_cum(:,jjb:jje,:)
	536	enddo
	537
	538	c=====================================================================
	539	c Transport méridien
	540	c=====================================================================
	541
	542	c cumul zonal des masses des mailles
	543	c ----------------------------------
[985]	544	jjb=jj_begin
	545	jje=jj_end
	546	if (pole_sud) jje=jj_end-1
	547
	548	zv(jjb:jje,:)=0.
	549	zmasse(jjb:jje,:)=0.
	550
	551	call Register_Hallo(masse_cum,ip1jmp1,llm,1,1,1,1,Req)
	552	do iQ=1,nQ
	553	call Register_Hallo(Q_cum(1,1,1,iQ),ip1jmp1,llm,0,1,1,0,Req)
	554	enddo
	555
	556	call SendRequest(Req)
	557	call WaitRequest(Req)
	558
	559	call massbar_p(masse_cum,massebx,masseby)
[630]	560
	561	jjb=jj_begin
	562	jje=jj_end
	563	if (pole_sud) jje=jj_end-1
	564
	565	do l=1,llm
	566	do j=jjb,jje
	567	do i=1,iim
	568	zmasse(j,l)=zmasse(j,l)+masseby(i,j,l)
	569	zv(j,l)=zv(j,l)+flux_v_cum(i,j,l)
	570	enddo
	571	zfactv(j,l)=cv(1,j)/zmasse(j,l)
	572	enddo
	573	enddo
	574
	575	c print*,'3OK'
	576	c --------------------------------------------------------------
	577	c calcul de la moyenne zonale du transport :
	578	c ------------------------------------------
	579	c
	580	c --
	581	c TOT : la circulation totale [ vq ]
	582	c
	583	c - -
	584	c MMC : mean meridional circulation [ v ] [ q ]
	585	c
	586	c ---- -- - -
	587	c TRS : transitoires [ v'q'] = [ vq ] - [ v q ]
	588	c
	589	c - * - * - - - -
	590	c STT : stationaires [ v q ] = [ v q ] - [ v ] [ q ]
	591	c
	592	c - -
	593	c on utilise aussi l'intermediaire TMP : [ v q ]
	594	c
	595	c la variable zfactv transforme un transport meridien cumule
	596	c en kg/s * unte-du-champ-transporte en m/s * unite-du-champ-transporte
	597	c
	598	c --------------------------------------------------------------
	599
	600
	601	c ----------------------------------------
	602	c Transport dans le plan latitude-altitude
	603	c ----------------------------------------
	604
[985]	605	jjb=jj_begin
	606	jje=jj_end
	607	if (pole_sud) jje=jj_end-1
	608
[630]	609	zvQ=0.
	610	psiQ=0.
	611	do iQ=1,nQ
	612	zvQtmp=0.
	613	do l=1,llm
	614	do j=jjb,jje
	615	c print*,'j,l,iQ=',j,l,iQ
	616	c Calcul des moyennes zonales du transort total et de zvQtmp
	617	do i=1,iim
	618	zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)
	619	s +flux_vQ_cum(i,j,l,iQ)
	620	zqy= 0.5(Q_cum(i,j,l,iQ)masse_cum(i,j,l)+
	621	s Q_cum(i,j+1,l,iQ)*masse_cum(i,j+1,l))
	622	zvQtmp(j,l)=zvQtmp(j,l)+flux_v_cum(i,j,l)*zqy
	623	s /(0.5*(masse_cum(i,j,l)+masse_cum(i,j+1,l)))
	624	zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)+zqy
	625	enddo
	626	c print*,'aOK'
	627	c Decomposition
	628	zvQ(j,l,iave,iQ)=zvQ(j,l,iave,iQ)/zmasse(j,l)
	629	zvQ(j,l,itot,iQ)=zvQ(j,l,itot,iQ)*zfactv(j,l)
	630	zvQtmp(j,l)=zvQtmp(j,l)*zfactv(j,l)
	631	zvQ(j,l,immc,iQ)=zv(j,l)zvQ(j,l,iave,iQ)zfactv(j,l)
	632	zvQ(j,l,itrs,iQ)=zvQ(j,l,itot,iQ)-zvQtmp(j,l)
	633	zvQ(j,l,istn,iQ)=zvQtmp(j,l)-zvQ(j,l,immc,iQ)
	634	enddo
	635	enddo
	636	c fonction de courant meridienne pour la quantite Q
	637	do l=llm,1,-1
	638	do j=jjb,jje
	639	psiQ(j,l,iQ)=psiQ(j,l+1,iQ)+zvQ(j,l,itot,iQ)
	640	enddo
	641	enddo
	642	enddo
	643
	644	c fonction de courant pour la circulation meridienne moyenne
[985]	645	psi(jjb:jje,:)=0.
[630]	646	do l=llm,1,-1
	647	do j=jjb,jje
	648	psi(j,l)=psi(j,l+1)+zv(j,l)
	649	zv(j,l)=zv(j,l)*zfactv(j,l)
	650	enddo
	651	enddo
	652
	653	c print*,'4OK'
	654	c sorties proprement dites
	655	if (i_sortie.eq.1) then
	656	jjb=jj_begin
	657	jje=jj_end
	658	jjn=jj_nb
	659	if (pole_sud) jje=jj_end-1
	660	if (pole_sud) jjn=jj_nb-1
	661
	662	do iQ=1,nQ
	663	do itr=1,ntr
	664	call histwrite(fileid,znom(itr,iQ),itau,
	665	s zvQ(jjb:jje,:,itr,iQ)
	666	s ,jjn*llm,ndex3d)
	667	enddo
	668	call histwrite(fileid,'psi'//nom(iQ),
	669	s itau,psiQ(jjb:jje,1:llm,iQ)
	670	s ,jjn*llm,ndex3d)
	671	enddo
	672
[985]	673	call histwrite(fileid,'masse',itau,zmasse(jjb:jje,1:llm)
[630]	674	s ,jjn*llm,ndex3d)
[985]	675	call histwrite(fileid,'v',itau,zv(jjb:jje,1:llm)
[630]	676	s ,jjn*llm,ndex3d)
[985]	677	psi(jjb:jje,:)=psi(jjb:jje,:)*1.e-9
[630]	678	call histwrite(fileid,'psi',itau,psi(jjb:jje,1:llm),
	679	s jjn*llm,ndex3d)
	680
	681	endif
	682
	683
	684	c -----------------
	685	c Moyenne verticale
	686	c -----------------
	687
[985]	688	zamasse(jjb:jje)=0.
[630]	689	do l=1,llm
	690	zamasse(jjb:jje)=zamasse(jjb:jje)+zmasse(jjb:jje,l)
	691	enddo
[985]	692
	693	zavQ(jjb:jje,:,:)=0.
[630]	694	do iQ=1,nQ
	695	do itr=2,ntr
	696	do l=1,llm
	697	zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ)
	698	s +zvQ(jjb:jje,l,itr,iQ)
	699	s *zmasse(jjb:jje,l)
	700	enddo
	701	zavQ(jjb:jje,itr,iQ)=zavQ(jjb:jje,itr,iQ)/zamasse(jjb:jje)
	702	call histwrite(fileid,'a'//znom(itr,iQ),itau,
	703	s zavQ(jjb:jje,itr,iQ),jjn*llm,ndex3d)
	704	enddo
	705	enddo
	706
	707	c on doit pouvoir tracer systematiquement la fonction de courant.
	708
	709	c=====================================================================
	710	c/////////////////////////////////////////////////////////////////////
	711	icum=0 !///////////////////////////////////////
	712	endif ! icum.eq.ncum !///////////////////////////////////////
	713	c/////////////////////////////////////////////////////////////////////
	714	c=====================================================================
	715
	716	return
	717	end

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

Download in other formats:

Original Format