Context Navigation

bilan_dyn_p.F @ 1080

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

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

Download in other formats:

Original Format