Context Navigation

bilan_dyn_p.F @ 1114

Last change on this file since 1114 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

Line
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
166	! definition du domaine d'ecriture pour le rebuild
167
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
179	c=====================================================================
180	c Initialisation
181	c=====================================================================
182	ndex3d=0
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
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)
262
263	call histbeg(trim(infile),
264	. 1, rlong(jjb:jje), jjn, rlatg(jjb:jje),
265	. 1, 1, 1, jjn,
266	. tau0, zjulian, dt_cum, thoriid, fileid,
267	. bilan_dyn_domain_id)
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,:)
372	Q(:,jjb:jje,:,iovap)=trac(:,jjb:jje,:,1)
373	Q(:,jjb:jje,:,iun)=1.
374
375
376	c=====================================================================
377	c Cumul
378	c=====================================================================
379	c
380	if(icum.EQ.0) then
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.
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
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
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
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
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
481	jjb=jj_begin
482	jje=jj_end
483
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
513	jjb=jj_begin
514	jje=jj_end
515
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
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 ----------------------------------
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)
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
603	jjb=jj_begin
604	jje=jj_end
605	if (pole_sud) jje=jj_end-1
606
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
643	psi(jjb:jje,:)=0.
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
671	call histwrite(fileid,'masse',itau,zmasse(jjb:jje,1:llm)
672	s ,jjn*llm,ndex3d)
673	call histwrite(fileid,'v',itau,zv(jjb:jje,1:llm)
674	s ,jjn*llm,ndex3d)
675	psi(jjb:jje,:)=psi(jjb:jje,:)*1.e-9
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
686	zamasse(jjb:jje)=0.
687	do l=1,llm
688	zamasse(jjb:jje)=zamasse(jjb:jje)+zmasse(jjb:jje,l)
689	enddo
690
691	zavQ(jjb:jje,:,:)=0.
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