Context Navigation

bilan_dyn_p.F @ 5434

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

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

Download in other formats:

Original Format